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1.  Introduction 


Data  assimilation  is  used  in  numerical  weather  prediction  (NWP)  to  incorporate 
observations  into  the  NWP  model  to  provide  a  more  accurate  estimate  of  the 
meteorological  conditions.  There  are  a  wide  variety  of  data  assimilation 
techniques  used  in  NWP,  each  with  its  own  advantages  and  disadvantages. 
Recently,  a  multitude  of  combination  techniques  that  aim  to  leverage  the  strengths 
of  2  techniques  while  minimizing  the  weaknesses  of  the  individual  techniques 
have  emerged.  In  this  study,  we  propose  such  a  combination  data  assimilation 
technique  and  demonstrate  initial  results  of  applying  the  combination  technique  to 
a  single  case  day. 

Data  assimilation  techniques  used  in  NWP  include  3-D  variational  (3DVAR;  e.g., 
Hu  et  al.  2006;  Xiao  and  Sun  2007),  4-D  variational  (4DVAR;  e.g.,  Mahfouf  et  al. 
2005;  Lopez  201 1),  Ensemble  Kalman  filter  (EnKF;  e.g.,  Snyder  and  Zhang  2003; 
Dowell  et  al.  2011),  and  Newtonian  relaxation  (nudging;  e.g.,  Schroeder  et  al. 
2006;  Reen  and  Stauffer  2010). 

In  theory,  3DVAR  could  estimate  an  optimal  analysis  for  the  provided 
background  field  and  observations.  However,  since  creating  an  optimal  analysis 
requires  perfect  knowledge  of  the  background-error  covariance  and  the 
observation-error  covariance,  the  actual  analysis  will  not  be  optimal  in  practice. 
Background-error  covariance  is  often  estimated  by  using  previous  NWP 
simulations  or  multiple  concurrent  NWP  simulations,  and  thus  3DVAR  is 
challenging  to  apply  for  on-demand  forecasts  where  no  previous  NWP 
simulations  necessarily  exist  for  the  desired  domain  and  computational  resources 
do  not  allow  for  multiple  concurrent  NWP  simulations.  Without  the  use  of 
multiple  concurrent  NWP  simulations  to  calculate  the  background  error 
covariance,  the  background  error  covariance  will  lack  flow  dependence.  3DVAR 
also  has  the  disadvantage  that  it  is  an  intermittent  technique,  meaning  that  an 
analysis  is  created  and  applied  at  a  discrete  time.  Direct  insertion  of  the  analysis  at 
a  single  time  can  result  in  a  solution  with  noise  as  the  model  attempts  to  adjust  to 
a  solution  that  is  dynamically  consistent  with  the  model’s  equations.  One  way  to 
mitigate  this  is  to  apply  the  increment  calculated  using  a  3DVAR  analysis  over  a 
time  period  using  Incremental  Analysis  Updating  as  in  Brewster  and  Stratman 
(2016)  or  to  analysis-nudge  toward  the  analysis.  A  weakness  of  all  intermittent 
data  assimilation  techniques  regards  the  application  of  observations  made  at  a 
time  other  than  the  time  for  which  the  data  assimilation  technique  is  creating  an 
analysis.  The  effect  of  these  observations  cannot  be  dependent  on  the  difference 
between  the  observations  and  the  model  state  at  the  time  of  the  observations.  First 
guess  at  appropriate  time  (e.g.,  Lorenc  and  Rawlins  2005)  attempts  to  mitigate 
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this  by  comparing  observations  to  background  fields  closer  in  time  to  the 
observations,  but  this  is  only  a  rough  approximation  of  a  full  continuous  method 
since  1)  the  comparison  is  either  against  an  analysis  or  output  from  a  previous 
model  integration  and  thus  is  not  the  current  model  integration’s  solution  at  the 
observation  time  (which  will  be  dependent  on  assimilation  of  earlier 
observations),  2)  the  comparison  is  limited  to  the  temporally  nearest  analysis  (or 
temporal  interpolation  between  the  nearest  analyses),  and  3)  the  effect  of  the 
observation  is  still  applied  at  the  analysis  time  rather  than  the  time  of  the 
observation. 

4DVAR  involves  determining  a  linearized  version  of  the  NWP  model  (e.g., 
tangent  linear  model  or  perturbation  forecast  model;  Lawless  et  al.  2003)  as  well 
as  its  adjoint  and  then  using  these  models  to  allow  for  forward  and  backward 
integration  to  determine  the  initial  condition  analysis  that  should  result  in  the  best 
forecast  over  the  time  period  4DVAR  is  applied.  4DVAR  is  basically  3DVAR 
with  a  numerical  forecast  model  as  strong  constraint  and  the  temporal  dimension 
added.  Techniques  that  account  for  the  temporal  dimension  (i.e.,  accounting  for 
the  temporal  variation  of  observations  in  creating  an  analysis  at  a  single  time)  but 
do  not  use  the  numerical  forecast  model  as  a  strong  constraint  are  sometimes 
referred  to  as  3.5DVAR  (e.g.,  Zhao  et  al.  2008).  While  the  addition  of  the 
temporal  dimension  in  4DVAR  should  allow  the  analysis  to  incorporate 
observation  information  more  effectively  over  a  time  range  than  3DVAR,  one 
challenge  involved  in  using  4DVAR  is  the  very  significant  computational 
demands  of  the  method.  Also,  while  4DVAR  could  theoretically  be  implemented 
by  using  the  full  nonlinear  forecast  model  for  forward  integration  and  its  adjoint 
(if  it  exists)  for  backward  integration,  because  of  the  computational  expense 
involved  in  running  the  full  nonlinear  forecast  model  and  the  difficulty  in  finding 
the  adjoint,  in  practice  this  procedure  is  often  modified.  Instead,  the  full  nonlinear 
forecast  model  is  integrated  forward,  and  this  prediction  is  compared  to 
observations.  This  prediction  is  used  with  a  linearized  version  of  the  model 
(tangent  linear  model)  for  forward  integration  and  the  adjoint  of  the  linearized 
model  for  backward  integration  in  an  iterative  manner  to  find  an  updated  version 
of  the  initial  conditions.  This  iteration  is  referred  to  as  the  inner  loop. 

In  addition  to  being  linearized,  the  forward  model  may  be  coarser  and  more 
simplified  than  the  full  nonlinear  forecast  model.  The  proposed  initial  conditions 
arising  from  the  inner  loop  can  then  be  used  to  create  another  solution  using  the 
full  nonlinear  forecast  model,  which  then  can  have  the  inner  loop  applied  to  it. 
The  repeated  use  of  the  full  nonlinear  forecast  model  is  referred  to  as  the  outer 
loop  and  is  expected  to  decrease  errors  resulting  from  the  linearization  in  the  inner 
loop.  However,  computational  requirements  limit  the  number  of  integrations  of 
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the  outer  loop,  and  thus  the  linearization  and  simplification  applied  in  the  inner 
loop  can  hinder  the  capability  of  4DVAR  to  fully  reach  its  potential. 

EnKF  relies  on  an  ensemble  of  NWP  model  runs  to  diagnose  the  background- 
error  covariance.  This  allows  observations  to  be  assimilated  using  flow-dependent 
background-error  covariances.  EnKF  is  expensive  computationally  since  by 
definition  multiple  NWP  simulations  must  be  carried  out.  Additionally,  since  a 
very  large  number  of  ensemble  members  is  needed  to  accurately  specify  the 
background-error  covariances,  techniques  (such  as  localization)  must  be  used  to 
account  for  the  limited  number  of  ensemble  members  that  can  be  used  in  practice. 

Newtonian  relaxation  (also  known  as  nudging)  involves  adding  nonphysical  terms 
to  the  model’s  tendency  terms  to  gradually  nudge  the  model  toward  observations 
or  an  analysis.  Nudging  has  the  advantages  that  it  is  a  continuous  data 
assimilation  method  and  is  not  computationally  intensive.  Additionally,  since  the 
parameters  used  to  control  the  strength  of  nudging  and  the  spatiotemporal 
influence  of  observations  are  generally  set  based  on  past  studies,  it  does  not 
require  additional  NWP  model  integrations  to  specify  how  the  effect  of 
observations  are  spread  for  a  specific  domain  or  case.  However,  if  resources 
allow,  experiments  to  improve  these  parameters  could  be  beneficial.  Nudging  has 
the  disadvantages  that  it  cannot  directly  assimilate  variables  that  are  not 
prognostic  variables  and  that  it  does  not  fully  account  for  case-specific 
background-error  covariances.  A  more  detailed  description  of  nudging  is  provided 
in  Section  2.2. 

The  variational  version  of  the  Local  Analysis  and  Prediction  System  (vLAPS; 
Jiang  et  al.  2015)  uses  a  4-D,  multiscale  variational  scheme  (3D VAR)  without  an 
actual  model  constraint  for  pressure,  temperature,  wind,  and  humidity,  and  uses  a 
modified  Barnes  scheme  with  some  1-D  variational  analysis  components  (Albers 
et  al.  1996)  for  other  variables.  Prior  to  the  introduction  of  vLAPS,  traditional 
LAPS  used  the  latter  approach  for  all  variables.  Compared  to  other  variational 
techniques,  vLAPS  advantages  include  the  multiscale  nature  of  the  technique  (in 
space  and  time)  and  its  computational  efficiency.  The  multiscale  nature  of  the 
technique  improves  the  analysis,  especially  in  cases  where  the  background  error 
covariance  is  not  well  specified;  this  is  important  because  it  is  difficult  to 
accurately  determine  the  background-error  covariance  (Xie  et  al.  2011).  A  more 
detailed  description  of  vLAPS  is  provided  in  Section  2.3. 

Various  combinations  of  data  assimilation  techniques  have  been  implemented. 
One  example  is  the  operational  Weather  Research  and  forecasting  (WRL)  model- 
based  9-km  Rapid  Refresh  model,  which  uses  a  technique  referred  to  as  hybrid 
ensemble-variational  (Benjamin  et  al.  2016).  The  method  uses  3DVAR  but 
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specifies  the  background-error  covariance  as  a  combination  of  static  background- 
error  covariance  and  background-error  covariance  determined  from  an 
80-member  ensemble  of  a  coarser  global  model. 

There  have  been  several  combination  data  assimilation  techniques  that  include 
nudging  as  one  of  the  techniques.  Lei  et  al.  (2012a,  2012b,  2012c)  developed  a 
hybrid  nudging-ensemble  Kalman  filter  (HNEnKF)  first  in  simplified  models  and 
then  applied  it  to  the  Applied  Research  version  of  the  WRF  (Skamarock  et  al. 
2008).  This  method  used  an  ensemble  of  WRF  simulations  to  determine  both  the 
strength  with  which  observation  nudging  was  applied  and  the  horizontal  and 
spatial  spreading  of  the  nudging  for  each  observation.  In  addition,  information 
from  the  ensemble  was  used  to  allow  innovations  in  one  variable  to  affect  another 
variable  (e.g.,  an  error  in  the  u-wind  component  could  result  in  nudging  of  the  v- 
wind  component).  This  was  applied  in  a  series  of  3-h  forecasts  to  produce 
analyses  to  force  atmospheric  transport  and  dispersion  modeling  for  a  single  case. 
The  3-h  HNEnKF  forecasts  improved  over  EnKF  in  terms  of  wind  direction  and 
over  observation  nudging  in  terms  of  temperature  and  relative  humidity  while  the 
performance  in  the  other  fields  is  less  conclusive.  Lei  et  al.  (2012c)  also  show  that 
the  noise  levels  in  the  HNEnKF  and  4-D  data  assimilation  approaches  are  much 
lower  than  in  the  EnKF  simulation. 

The  National  Center  for  Atmospheric  Research  (NCAR)  developed  a  hybrid 
nudging-ensemble  system  they  refer  to  as  4D  Relaxation  Ensemble  Kalman  Filter 
(Liu  et  al.  2015).  It  uses  an  ensemble  to  determine  the  nudging  coefficients  used 
in  observation  nudging. 

Shaw  et  al.  (2008)  used  3DVAR  to  create  an  analysis  to  start  WRF  and  then 
performed  observation  nudging  during  the  next  3  h.  They  briefly  showed  results 
from  a  single  case  suggesting  the  combination  of  3DVAR  and  nudging  performed 
better  than  either  method  separately. 

Lei  and  Hacker  (2015)  investigated  observation  nudging,  EnKF,  and  HNEnKF  in 
a  Lorenz  model.  Lorenz  models  involve  a  much  simpler  set  of  equations  than 
NWP  models,  such  as  WRF,  but  include  some  behaviors  reminiscent  of  the 
atmosphere;  they  can  be  used  to  test  assimilation  techniques  before  application  to 
the  much  more  complex  environment  of  a  NWP  model.  They  find  that  nudging 
performs  better  with  increasing  model  error,  but  the  increased  model  error 
degrades  the  covariances  used  in  the  hybrid  scheme.  They  conclude  that  in  the 
Lorenz  model  they  investigated,  the  hybrid  scheme  cannot  result  in  errors  that  are 
simultaneously  lower  than  both  EnKF  and  observation  nudging.  In  other  words,  in 
the  Lorenz  model  for  any  individual  case  the  hybrid  scheme  cannot  outperform 
using  EnKF  alone  while  also  outperforming  using  observation  nudging  alone. 


Approved  for  public  release;  distribution  is  unlimited. 

4 


However,  it  is  not  clear  if  this  conclusion  holds  for  much  more  complicated 
models  such  as  WRF.  Additionally,  even  if  that  hybrid  method  (HNEnKF)  never 
performs  better  than  both  component  data  assimilation  methodologies,  the 
methodology  can  still  be  valuable  if  the  hybrid  method  normally  performs  better 
than  one  of  the  components  and  which  component  performs  best  varies  by  case; 
this  can  allow  the  hybrid  method  to  avoid  the  worst  performances  of  both 
methods,  which  might  be  exactly  the  cases  where  it  is  most  important  that 
forecasts  be  improved.  Additionally,  this  result  was  specifically  for  an  EnKF 
observation  nudging  hybrid  and  thus  is  not  directly  applicable  to  the  combination 
assimilation  methodologies  investigated  in  this  study. 

This  study  demonstrates  the  combination  of  the  vLAPS  and  nudging  data 
assimilation  techniques  for  a  single  case  day  with  strong  convection.  Section  2 
describes  the  Advanced  Research  version  of  the  WRF  model  (WRF-ARW)  and  its 
configuration  here,  as  well  as  describing  the  vLAPS  and  nudging  data 
assimilation  techniques.  A  description  of  the  case  is  given  in  Section  3,  Section  4 
details  the  methodology  used  to  carry  out  the  hybrid  data  assimilation,  and 
Section  5  describes  the  experimental  design.  Section  6  describes  preliminary 
results,  and  Section  7  contains  summary  and  a  discussion  of  the  future  work. 

2.  Model  Description  and  Configuration 


2.1  WRF-ARW 

The  Advanced  Research  version  of  the  Weather  Research  and  Forecasting  model 
(WRF-ARW)  V3.6.1  (Skamarock  et  al.  2008)  is  applied  with  56  vertical  layers 
for  an  801  x  801  1-km  horizontal  grid  spacing  domain  centered  over  Oklahoma 

(Fig-  1)- 
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Fig.  1  Area  covered  by  the  WRF  1-km  domain 

WRF-ARW  is  configured  here  to  use  the  Mellor-Yamada-Janjic  (MYJ)  scheme 
(Janjic  2001)  to  parameterize  the  atmospheric  boundary  layer  (ABL).  It  predicts 
turbulent  kinetic  energy  (TKE)  and  is  a  Mellor-Yamada  Level  2.5  turbulence 
closure  model.  As  in  Lee  et  al.  (2012)  and  Reen  et  al.  (2014),  the  background 
TKE  is  decreased  to  better  simulate  conditions  with  low  TKE,  and  the  ABL  depth 
diagnosis  is  altered. 

The  Thompson  microphysics  parameterization  (Thompson  et  al.  2008)  is  used, 
but  no  cumulus  parameterization  is  employed  because  of  the  high-resolution 
nature  of  the  single  domain.  For  radiation,  the  Rapid  Radiative  Transfer  Model 
(Mlawer  et  al.  1997)  is  used  for  longwave  and  the  Dudhia  scheme  (Dudhia  1989) 
for  shortwave.  Land  surface  processes  are  represented  by  the  Noah  land  surface 
model  (Chen  and  Dudhia  2001). 
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2.2  Nudging  Data  Assimilation 


Nudging  (e.g.,  Stauffer  and  Seaman  1990;  Stauffer  et  al.  1991),  also  known  as 
Newtonian  relaxation,  involves  nudging  the  model  toward  observations 
(observation  nudging)  or  analyses  (analysis  nudging)  by  adding  a  term  to  the 
tendency  equations  whose  magnitude  depends  on  the  magnitude  of  the  difference 
between  the  model  and  either  the  observation  or  the  analysis.  WRF  has  the  ability 
to  nudge  toward  potential  temperature,  water  vapor  mixing  ratio,  and  winds. 

Observation  nudging  (e.g.,  Stauffer  and  Seaman  1994;  Reen  2016)  calculates  the 
difference  between  an  observation  and  the  current  model  solution  at  that  point 
(the  difference  is  called  the  innovation)  and  then  uses  the  innovation  to  calculate 
an  additional  term  for  the  tendency  equation.  Observation  nudging  can  be 
described  by  the  following  equation  and  its  description  taken  directly  from  Reen 
(2016): 

N 

d  ^wg%x>y>z>tho(i)- 

-^-{pc,y,z,t)  =  Fq(x,y,z,t)+  juGq^ - * - ,  (1) 

YjWq(i,x,y,z,t) 

i= 1 

where  q  is  the  quantity  being  nudged  (e.g.,  water  vapor  mixing  ratio),  y  is  the  dry 
hydrostatic  pressure,  Fq  represents  the  physical  tendency  terms  of  q  (for  water 
vapor  this  includes  advection,  diffusion,  conversion  from  water  vapor  to  cloud 
water,  etc.),  Gq  is  the  nudging  strength  for  q,  N  is  the  total  number  of 
observations,  i  is  the  index  to  the  current  observation,  Wq  is  the  spatiotemporal 
weighting  function  based  on  the  temporal  and  spatial  separation  between  the 
observation  and  the  current  model  location,  q0  is  the  observed  value  of  q,  and 
qm(xuyi,zi,t )  is  the  model  value  of  q  interpolated  to  the  observation  location.  The 
quantity  q0-qm  is  the  innovation;  the  innovation  associated  with  a  given 
observation  evolves  with  time  (both  before  and  after  the  time  of  the  observation) 
as  the  model  value  (qm)  evolves.  Thus,  as  the  model  value  approaches  the 
observed  value,  the  nudging  tendency  term  decreases. 

To  apply  observation  nudging  effectively,  the  observation  must  be  at  least 
somewhat  representative  of  scales  represented  by  the  size  of  the  grid  cells  of  the 
model  integration  it  is  being  applied  to.  This  should  not  be  problematic  in  this 
study  because  observation  nudging  has  been  successfully  applied  in  past  studies 
to  horizontal  grid  spacings  much  coarser  than  the  1-km  horizontal  grid  spacing 
used  here.  The  strength  of  the  observation  nudging  used  here  ( Gq )  is  6  x  10  4  s  ', 
which  means  that  the  e-folding  time  for  the  error  is  approximately  28  min.  Thus, 
assuming  that  the  only  nonzero  tendency  term  for  the  observed  value  is  the 
observation  nudging  term,  the  difference  between  the  model  and  the  observation 
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after  approximately  28  min  would  be  the  previous  difference  divided  by  e 
(i.e.,  ~37%  of  the  previous  difference). 

The  innovation  for  an  observation  is  not  applied  at  a  single  time  but  rather  over  a 
time  period  since  the  errors  at  the  observation  location  are  correlated  with  errors 
nearby  in  time  (the  time  weighting  is  included  in  Wq  in  Eq.  1).  This  also  allows 
the  model  physical  tendency  terms  to  dominate,  and  thus  the  observations  can  be 
assimilated  while  maintaining  a  physically  consistent  solution  with  a  minimum  of 
“noise”.  Here,  full  weighting  of  above-surface  observations  is  applied  for  a  1.5-h 
time  period  centered  on  the  observation  time,  and  the  weighting  decreases  to  zero 
in  the  0.75-h  time  period  on  either  side  of  this  time  period.  For  surface 
observations  the  temporal  window  is  two-thirds  as  large  (full  strength  for  1.0  h 
and  decreasing  to  zero  in  the  0.5  h  on  either  side  of  this  time  period)  to  account 
for  the  more  rapid  atmospheric  changes  near  the  surface.  The  innovation  is 
recalculated  frequently  throughout  the  time  period  over  which  the  innovation  is 
applied  so  that  nudging  does  not  “overshoot”  and  start  nudging  the  model  away 
from  the  observation  (here  it  is  recalculated  every  6  s,  which  is  every  other  time 
step). 

The  innovation  is  also  applied  over  a  horizontal  spatial  extent  since  the  errors  at 
the  observation  location  are  correlated  with  errors  nearby  spatially  (the  spatial 
weighting  is  included  in  Wq  in  Eq.  1).  In  general,  the  innovations  are  spread 
horizontally  with  isotropic  weighting  functions,  with  weight  decreasing  with 
distance  from  the  observations.  Here,  a  30-km  radius  of  influence  is  used  at  the 
surface,  60  km  just  above  the  surface,  and  the  radius  of  influence  increases  with 
decreasing  pressure  up  to  a  maximum  of  120  km  at  500  hPa.  Innovations  from 
surface  observations  are  spread  along  the  surface,  but  the  strength  at  which  the 
innovation  is  applied  is  decreased  for  areas  whose  model  terrain  height  differs 
from  that  at  the  location  of  the  observation.  This  is  used  to  limit  the  spreading  of 
surface  observations  in  a  valley  to  a  mountaintop  and  vice  versa. 

The  vertical  spreading  of  innovations  depends  on  the  observation  type  (the 
vertical  spreading  is  also  included  in  Wq  in  Eq.  1).  For  multilevel  observations 
(e.g.,  radiosondes),  the  innovation  profile  at  the  location  of  the  observation  is 
interpolated  vertically  to  the  model  vertical  levels  at  the  horizontal  location  that 
the  innovation  is  being  applied.  For  single-level  above-surface  observations  (e.g., 
aircraft  observations),  the  innovation  is  applied  between  50  hPa  below  and  50  hPa 
above  the  observation  with  the  strength  at  which  it  is  applied  decreasing  linearly 
with  pressure  difference  between  the  observation  and  the  vertical  location  at 
which  the  innovation  is  being  applied.  Additionally,  the  innovations  calculated 
above  the  top  of  the  model-diagnosed  ABL  are  not  applied  below  the  top  of  the 
model-diagnosed  ABL,  and  vice  versa.  For  surface  observations,  in  convective 
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conditions  the  innovation  is  applied  with  full  weight  throughout  the  ABL  and  then 
ramps  down  to  zero  in  the  next  50  m.  In  stable  conditions  the  innovation  is 
applied  with  full  weight  in  the  lowest  50  m  and  then  ramps  down  to  zero  in  the 
next  50  m. 

There  are  some  additional  details  regarding  the  observation  nudging 
configuration.  The  modification  of  observation  nudging  of  water  vapor  by  Reen  et 
al.  (2016)  is  applied  here  to  prevent  excessive  drying.  At  the  end  of  the 
observation  nudging  period,  we  use  a  1-h  ramp-down  period  to  gradually  reduce 
the  observation  nudging  tendency  term  to  minimize  any  noise  that  might  be 
created  if  the  observation  nudging  tendency  terms  were  suddenly  removed. 
During  this  ramp-down  period,  2  things  differ  from  observation  nudging  prior  to 
this  ramp-down  period.  First,  no  new  observations  (i.e.,  observations  taken  during 
the  ramp-down  period)  are  assimilated,  thus  nudging  is  toward  observations  taken 
prior  to  the  ramp-down  period  that  are  not  yet  too  old  to  be  outside  the  temporal 
window  of  that  observation.  Second,  a  linear  ramp-down  weighting  factor  is 
multiplied  by  the  other  weighting  factors  to  gradually  ramp  down  the  observation 
nudging  during  this  period.  An  unexpectedly  large  computational  expense  was 
associated  with  observation  nudging  in  these  experiments.  Further  investigation  is 
needed  to  determine  the  causes  of  this  expense.  It  may  be  that  the  large  number  of 
model  grid  cells  in  the  domain  (-640,000  per  level),  compared  to  most 
applications  of  observation  nudging,  exposed  an  area  where  improvements  are 
needed  in  the  implementation  of  observation  nudging  in  WRF  to  allow  it  to  better 
scale  to  larger  domains. 

Analysis  nudging  (e.g.,  Stauffer  and  Seaman  1994;  Gilliam  et  al.  2012)  calculates 
the  difference  between  an  analysis  and  the  model  and  then  uses  that  difference 
(innovation)  to  construct  an  additional  term  for  the  tendency  equations  to  nudge 
the  model  toward  the  analysis.  To  apply  analysis  nudging  effectively,  the  analysis 
must  be  representative  of  the  scales  resolved  by  the  model  domain  being  nudged 
toward.  If  the  analysis  only  resolves  scales  much  coarser  than  the  model  domain, 
then  analysis  nudging  will  work  to  remove  the  finer-scale  structures  the  model 
simulates;  there  is  an  alternative  form  of  analysis  nudging  known  as  spectral 
nudging  that  nudges  only  the  large-scale  structures  (e.g.,  Vincent  and  Hahmann 
2015).  Analysis  nudging  is  usually  applied  to  relatively  coarse  domains  (e.g., 
Stauffer  and  Seaman  1994;  Rogers  et  al.  2013;  Exposito  et  al.  2015)  since 
analyses  do  not  normally  have  sufficient  coverage  of  high-resolution  observations 
to  make  the  analysis  representative  of  finer  domains.  However,  in  this  study  the 
use  of  high-resolution  data  sources  in  the  creation  of  analyses  (as  is  described 
later)  allow  analysis  nudging  to  be  applied  on  the  1-km  WRF  domain.  The 
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strength  of  the  analysis  nudging  used  here  is  3  x  10  4  s  ',  which  means  that  the  e- 
folding  time  for  the  error  is  approximately  56  min. 

Since  analysis  nudging  is  toward  an  analysis  with  values  at  each  model  grid  point 
at  specific  time  intervals,  the  temporal  and  spatial  spreading  applied  differs  from 
that  used  in  observation  nudging.  The  tendency  terms  added  due  to  analysis 
nudging  are  applied  at  all  model  time  steps  during  the  data  assimilation  period. 
Analyses  are  temporally  interpolated  to  the  current  time  to  determine  the  analysis 
against  which  the  innovation  will  be  calculated.  The  analysis  nudging  innovation 
is  calculated  here  at  each  time  step.  Spatially,  each  model  grid  point  is  generally 
nudged  based  on  the  coincident  analysis  grid  point,  although  this  differs  near  the 
surface.  WRF  analysis  nudging  also  has  the  capability  to  limit  the  strength  of 
analysis  nudging  in  locations  where  limited  nearby  observations  were  included  in 
the  analysis;  in  this  study,  because  of  the  nature  of  the  analyses,  they  were 
assumed  to  be  equally  valid  at  all  model  grid  points. 

Here,  we  nudge  both  toward  surface  analyses  (surface  analysis  nudging)  and 
analyses  above  the  surface  (3-D  analysis  nudging).  When  both  surface  and  3-D 
analysis  nudging  are  activated,  1)  the  surface  analyses  innovations  are  applied 
with  full  strength  throughout  the  ABL  and  with  90%  of  full  weight  in  the  layer 
immediately  above  the  ABL  and  2)  the  above-surface  analyses  are  applied  at  full 
strength  at  all  levels  above  the  top  of  the  ABL  except  at  10%  of  full  weight  in  the 
layer  immediately  above  the  ABL. 

While  surface  analysis  nudging  of  potential  temperature  is  straightforward, 
additional  details  are  involved  with  wind  and  moisture.  For  potential  temperature, 
the  model’s  lowest  prognostic  level  temperature  is  compared  to  the  analysis,  and 
the  resulting  innovation  applied  throughout  the  ABL.  Similarity  theory  is  used  to 
adjust  the  10-m  wind  analysis  to  the  lowest  model  layer  (e.g.,  Stauffer  et  al.  1991) 
if  the  top  of  the  ABL  is  above  the  lowest  model  layer  and  the  lowest  model  layer 
where  wind  is  predicted  is  above  10-m  AGL.  The  adjusted  surface  wind  u  and  v 
components  are  then  compared  to  the  model’s  lowest  prognostic  level  values  of  u 
and  v,  and  the  resulting  innovations  applied  throughout  the  ABL.  The  surface 
water  vapor  observation  is  compared  to  the  model  water  vapor  at  each  model 
prognostic  layer  in  the  ABL  (and  one  layer  above),  and  the  resultant  innovation 
applied  to  that  same  level. 

At  the  end  of  the  analysis  nudging  period,  we  use  a  15 -min  time  period  to  ramp 
down  the  influence  of  the  analysis  nudging.  During  this  time  period,  the  model  is 
nudged  toward  the  analysis  valid  at  the  beginning  of  this  time  period,  and  the 
weighting  decreases  linearly  with  time  to  zero.  Although  the  standard  version  of 
WRF  V3.6.1  does  not  nudge  toward  the  analysis  valid  at  the  beginning  of  the 
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ramp-down  time  period  because  of  a  bug,  the  bug  was  fixed  in  the  version  of 
WRF  V3.6.1  used  for  this  study  (the  bug  fix  was  also  submitted  to  the  WRF 
developers  and  included  in  the  public  release  version  of  WRF  starting  in  WRF 
V3.8). 

2.3  vLAPS  Data  Assimilation 

LAPS  was  developed  by  the  National  Oceanic  and  Atmospheric  Administration 
(McGinley  et  al.  1991)  and  traditionally  used  a  modified  Barnes  analysis 
(Hiemstra  et  al.  2006)  with  some  1-D  variational  components  (Albers  et  al.  1996). 
The  vLAPS  retains  the  traditional  techniques  for  some  variables  but  moves  to  a  3- 
D,  multiscale  variational  scheme  (3DVAR)  for  pressure,  temperature,  winds,  and 
humidity  (Jiang  et  al.  2015).  Some  of  the  vLAPS  description  in  this  section  is 
from  Jiang  et  al.  (2015). 

The  3DVAR  scheme  used  in  vLAPS  is  based  on  the  Space  and  Time  Multiscale 
Analysis  System  (Xie  et  al.  2011).  Even  though  it  does  not  have  a  forecast  model 
as  a  constraint  like  a  4DVAR,  it  indeed  considers  temporal  information  with  a 
half-hour  analysis  time  window.  The  multiscale  approach  pioneered  in 
computational  mathematics  applications  (see  Briggs  1987)  was  adopted  by  Xie  et 
al.  (2011)  for  meteorological  and  oceanographic  data  assimilation  and  improves 
accuracy  and  computational  efficiency  compared  to  nonmultiscale  variational 
techniques.  vLAPS  applies  the  multiscale  technique  spatially  and  temporally. 

Spatially,  vLAPS  completes  a  3DVAR  analysis  on  a  coarse  grid,  and  the  output  of 
that  analysis  is  used  as  the  background  for  a  3DVAR  analysis  at  a  finer  scale;  this 
process  continues  until  the  analysis  is  at  an  adequate  resolution.  For  each  grid,  a 
Laplacian  function  of  the  difference  between  analysis  and  background  is  added  to 
the  cost  function  to  mimic  the  background  term  of  a  3DVAR,  which  determines 
the  spreading  of  the  influence  of  an  observation.  The  coefficients  multiplying  the 
Laplacian  terms  control  the  influence  radius  of  observations.  The  use  of  a 
Laplacian  rather  than  a  recursive  filter  has  the  benefit  of  better  retaining  the 
gradients  of  the  background  to  fit  the  observations.  Since  vLAPS  reduces  the 
coefficients  of  the  Laplacian  terms  in  the  cost  function  when  moving  from  a 
coarser  grid  to  a  finer  one,  the  effective  radius  of  influence  of  an  observation 
decreases  as  the  horizontal  grid  spacing  of  the  analysis  decreases.  Performing  the 
analysis  at  multiple  spatial  scales  allows  the  influence  of  observations  in  data- 
sparse  areas  to  spread  broadly  in  coarser  grids  while  keeping  the  fine-scale 
features  resolvable  in  data-rich  areas  in  finer  grids.  This  also  makes  the  technique 
less  dependent  on  accurate  background  error  covariance.  Here,  the  coarsest  grid 
used  for  the  analysis  has  16-km  horizontal  grid  spacing  and  50-hPa  vertical  grid 
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spacing;  each  subsequent  analysis  decreases  the  horizontal  and  vertical  grid 
spacing  in  half  until  this  division  is  stopped  if  the  grid  space  exceeds  the  smallest 
spacing  specified  by  the  user.  For  this  case,  we  set  25  hPa  for  the  vertical.  This, 
the  finest  analysis  grid,  uses  2  km  horizontal  grid  spacing  and  25  vertical  grid 
spacing.  Five  analyses  were  completed,  including  the  coarsest  and  finest 
resolution  analyses. 

The  surface  analysis  is  completed  separately  from  the  analysis  for  layers  above 
the  surface,  and  the  surface  analysis  lacks  a  vertical  component. 

Temporally,  vLAPS  uses  multiple  time  frames  to  create  the  analysis  at  an 
individual  analysis  time;  this  allows  extraction  of  information  from  frequent 
observations  (e.g.,  from  radars)  related  to  rapid  changes.  A  time  filter  is  used  to 
weigh  the  contributions  of  observations  nearer  to  the  analysis  time  more  strongly 
than  observations  farther  in  the  past  or  future.  The  time  filter  used  is  a  Laplacian 
filter.  As  with  the  spatial  dimension,  the  analysis  is  repeated  multiple  times  with 
increasingly  short  temporal  windows  in  which  observations  can  affect  the 
analysis.  In  this  case,  this  results  in  the  longest  temporal  window  being  30  min, 
and  the  shortest  being  15  min,  both  centered  on  the  analysis  time.  Note  that  the 
spatial  and  temporal  refinement  of  the  analyses  are  taking  place  simultaneously 
(i.e.,  the  first  analysis  is  coarse  spatially  and  temporally,  the  next  analysis  is  finer 
spatially  and  temporally).  For  observations  that  are  not  taken  at  the  time  of  one  of 
the  first-guess  analyses,  the  first-guess  analyses  are  temporally  interpolated  to  the 
observation  to  determine  the  difference  between  the  observation  and  the  analyses. 

The  multiscale  approach  of  vLAPS  is  related  to  other  data  assimilation 
approaches.  For  example,  it  is  similar  to  a  4DVAR  approach,  but  instead  of  using 
the  full  NWP  model  as  a  constraint,  vLAPS  uses  a  time  filter  and  other  simplified 
constraints.  The  simplified  constraints  vLAPS  used  include  geostrophic  and 
hydrostatic  balances  as  weak  constraints  and  incompressibility  as  a  strong 
constraint.  It  is  also  similar  to  ensemble-variational  data  assimilation  (4DEnVar), 
wherein  ensemble  covariances  are  used  for  temporal  constraints  (Liu  et  al.  2008). 

The  vLAPS  variational  assimilation  differs  significantly  from  the  3DVAR  in  the 
widely  used  Gridpoint  Statistical  Interpolation  (GSI)  package.  In  addition  to  the 
multigrid  technique,  for  control  variables  vLAPS  uses  the  u-  and  v-wind 
components  (Xie  and  MacDonald  2011),  whereas  GSI  uses  stream  function  and 
velocity  potential.  Xie  et  al.  (2002  and  2011)  demonstrate  that  the  use  of  the  latter 
2  for  control  variables  is  undesirable;  this  is  especially  true  for  convective-scale 
resolutions.  There  are  several  other  advantages  of  vLAPS,  such  as  better  use  of 
physical  balances,  hot-start,  more  efficient  storage  of  the  background  covariances, 
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nonlinear  minimization  compared  to  GSI’s  linearized  minimization,  and  better 
representation  of  observation  errors. 

The  cloud  analysis  technique  applied  here  was  first  described  in  Albers  et  al. 
(1996).  Geostationary  Operational  Environmental  Satellite  (GOES)  data  (infrared 
+  visible),  radar,  surface  observations  (METARs),  and  model  first-guess  fields  are 
combined  in  the  so-called  “traditional”  or  sequential  cloud  analysis.  A 
preliminary  analysis  is  constructed  of  cloud  fraction  using  model  first-guess 
humidity  and  METAR  reports  that  are  spread  horizontally.  Next,  infrared  satellite 
data  are  added  in  a  series  of  cloud-clearing  and  cloud-building  steps,  with  similar 
steps  for  visible  satellite  data.  Finally,  radar  data  are  added  to  complete  the  3-D 
cloud  fraction  analysis.  This  cloud  field  is  postprocessed  into  3-D  hydrometeor 
fields  and  associated  vertical  velocity  to  feed  back  into  the  dynamical  fields 
(namely,  the  3-D  wind).  A  variational  cloud  analysis  is  under  development  for 
future  use. 

3.  Case  Description 

To  demonstrate  the  potential  utility  of  the  combined  data  assimilation  technique 
proposed  here,  we  investigate  a  single  case  day  (20  May  2013)  with  strong 
convection  in  the  southern  Great  Plains  region.  At  12  Coordinated  Universal  Time 
(UTC)  at  500  hPa  (Fig.  2),  an  upper-level  low  was  centered  over  the  northern 
Great  Plains  with  the  axis  of  the  trough  to  the  southwest.  This  resulted  in 
approximately  50-kt  (~26  ms'1)  flow  from  the  west-southwest  over  Oklahoma 
and  eastern  Kansas  at  500  hPa. 
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Fig.  2  Winds,  heights,  and  temperature  at  500  hPa  on  20  May  2013  at  12  UTC  from  the 
National  Center  for  Environmental  Prediction  daily  weather  maps 
(http://www.wpc.ncep.noaa.gov/dailywxmap).  Barbs  are  used  to  show  wind  (knots),  solid 
lines  show  isoheights  (dekameters  above  sea  level),  and  dashed  lines  show  isotherms  (°C). 

At  the  surface,  at  12  UTC  (Fig.  3)  a  cold  front  extended  south  from  the  northern 
Great  Plains  but  transitioned  to  a  stationary  front  over  Kansas  and  Oklahoma  and 
curved  more  westerly,  with  a  dry  line  extending  south  from  the  front  over  west 
Texas.  Observations  (Fig.  4a)  indicate  a  strong  moisture  contrast  near  the  front  in 
southern  Oklahoma  and  northern  Texas  at  12  UTC,  with  much  drier  conditions  on 
the  northwest  side  of  the  front.  Precipitation  at  12  UTC  (Fig.  4b)  is  limited  to  a 
broken  swath  oriented  west-southwest  to  east-northeast  from  far  northern  Texas 
through  northern  Oklahoma  and  southern  Kansas.  As  diurnal  heating  progresses, 
the  frontal  position  and  moisture  gradient  remain  similar  through  15  UTC 
(Fig.  4c)  and  18  UTC  (Fig.  4e)  with  the  weak  line  of  precipitation  becoming  more 
continuous  and  oriented  southwest  to  northeast  mostly  over  Kansas  (Fig.  4e-f). 

By  2100  UTC  a  dry  line  is  analyzed  south  from  near  the  front  in  central 
Oklahoma  (Fig.  4g)  with  a  line  of  strong  convection  in  place  through  central 
Oklahoma  and  southeastern  Kansas  (Fig.  4h).  At  0000  UTC  the  surface  analysis 
(Fig.  4i)  indicates  frontal  placement  is  similar,  but  now  analyzed  as  a  cold  front, 
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and  a  squall  line  is  analyzed  across  southeastern  Oklahoma,  which  is  consistent 
with  the  radar  observation  at  this  time  (Fig.  4j). 


Fig.  3  Surface  analyses  of  the  continental  United  States  from  the  Weather  Prediction 
Center  (http://www.wpc.ncep.noaa.gov/html/sfc2.shtml)  for  1200  UTC  20  May  2013 
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Fig.  4  Surface  analyses  (a,c,e,g,i)  of  the  southern  Great  Plains  from  the  Weather 
Prediction  Center  (http://www.wpc.ncep.noaa.gov/html/sfc2.shtml)  and  composite  base 
reflectivity  (b,d,f,h,j)  from  Iowa  Environmental  Mesonet  (http://mesonet.agron.iastate.edu/ 
GIS/apps/rview/warnings.phtml)  for  (a,b)  1200  UTC  20  May  2013,  (c,d)  1500  UTC  20  May 
2013,  (e,f)  1800  UTC  20  May  2013,  (g,h)  2100  UTC  20  May  2013,  and  (i,j)  0000  UTC  21  May 
2013 
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c) 


d) 


Fig.  4  Surface  analyses  (a,c,e,g,i)  of  the  southern  Great  Plains  from  the  Weather 
Prediction  Center  (http://www.wpc.ncep.noaa.gov/html/sfc2.shtml)  and  composite  base 
reflectivity  (b,d,f,h,j)  from  Iowa  Environmental  Mesonet  (http://mesonet.agron.iastate.edu/ 
GIS/apps  /rview/warnings.phtml)  for  (a,b)  1200  UTC  20  May  2013,  (c,d)  1500  UTC  20  May 
2013,  (e,f)  1800  UTC  20  May  2013,  (g,h)  2100  UTC  20  May  2013,  and  (i,j)  0000  UTC  21  May 
2013  (continued) 
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Fig.  4  Surface  analyses  (a,c,e,g,i)  of  the  southern  Great  Plains  from  the  Weather 
Prediction  Center  (http://www.wpc.ncep.noaa.gov/html/sfc2.shtml)  and  composite  base 
reflectivity  (b,d,f,h,j)  from  Iowa  Environmental  Mesonet  (http://mesonet.agron.iastate.edu/ 
GIS/apps  / rview/warnings.phtml)  for  (a,b)  1200  UTC  20  May  2013,  (c,d)  1500  UTC  20  May 
2013,  (e,f)  1800  UTC  20  May  2013,  (g,h)  2100  UTC  20  May  2013,  and  (i,j)  0000  UTC  21  May 
2013  (continued) 
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Fig.  4  Surface  analyses  (a,c,e,g,i)  of  the  southern  Great  Plains  from  the  Weather 
Prediction  Center  (http://www.wpc.ncep.noaa.gov/html/sfc2.shtml)  and  composite  base 
reflectivity  (b,d,f,h,j)  from  Iowa  Environmental  Mesonet  (http://mesonet.agron.iastate.edu/ 
GIS/apps  /rview/warnings.phtml)  for  (a,b)  1200  UTC  20  May  2013,  (c,d)  1500  UTC  20  May 
2013,  (e,f)  1800  UTC  20  May  2013,  (g,h)  2100  UTC  20  May  2013,  and  (i,j)  0000  UTC  21  May 
2013  (continued) 
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i) 


j) 


Fig.  4  Surface  analyses  (a,c,e,g,i)  of  the  southern  Great  Plains  from  the  Weather 
Prediction  Center  (http://www.wpc.ncep.noaa.gov/html/sfc2.shtml)  and  composite  base 
reflectivity  (b,d,f,h,j)  from  Iowa  Environmental  Mesonet  (http://mesonet.agron.iastate.edu/ 
GIS/apps  /rview/warnings.phtml)  for  (a,b)  1200  UTC  20  May  2013,  (c,d)  1500  UTC  20  May 
2013,  (e,f)  1800  UTC  20  May  2013,  (g,h)  2100  UTC  20  May  2013,  and  (i,j)  0000  UTC  21  May 
2013  (continued) 
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The  1200  UTC  radiosonde  at  Norman,  Oklahoma  (Fig.  5a),  indicates  a  thin 
saturated  layer  at  approximately  900  hPa  at  the  base  of  a  strong  capping  inversion 
and  then  a  deep  elevated  mixed  layer.  At  1800  UTC  (Fig.  5b),  daytime  heating 
along  with  the  formation  of  a  convective  boundary  layer  and  a  weakened  capping 
inversion  are  present;  the  convective  available  potential  energy  at  this  time  is 
3135  J  kg'1.  By  0000  UTC,  the  main  line  of  convection  is  southeast  of  Norman, 
but  an  isolated  cell  is  at  least  in  the  vicinity  of  Norman  (Fig.  4j),  and  there  is  a 
shallow  isothermal  layer  at  the  surface  topped  by  an  inversion  (Fig.  5c).  Also, 
note  the  backing  in  the  near-surface  wind  between  1800  and  0000  UTC  (compare 
Fig.  5c  to  Fig.  5b),  and  the  speed  shear  present  in  the  1800  and  0000  UTC 
radiosondes. 

A  total  of  42  tornadoes  were  reported  in  the  United  States  on  20  May  2013, 
including  18  in  Oklahoma,  11  in  Missouri,  5  in  Texas,  3  in  Kansas,  and  3  in 
Arkansas  (NCEI  2016).  An  EF5  touched  down  at  1956  UTC  in  Newcastle, 
Oklahoma,  and  traveled  through  Moore,  Oklahoma,  before  dissipating 
approximately  2035  UTC  (Atkins  et  al.  2014;  Burgess  et  al.  2014;  Kurdzo  et  al. 
2015). 


Fig.  5  Skew-T’s  of  radiosondes  from  Norman,  Oklahoma,  for  a)  12  UTC  20  May  2013,  b) 
18  UTC  20  May  2013,  and  c)  00  UTC  21  May  2013  (obtained  from  http://weather.uwyo.edu 
/upperair/sounding.html).  A  full  wind  barb  represents  10  kts  and  a  half  wind  barb 
represents  5  kts. 
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Fig.  5  Skew-T’s  of  radiosondes  from  Norman,  Oklahoma,  for  a)  12  UTC  20  May  2013,  b) 
18  UTC  20  May  2013,  and  c)  00  UTC  21  May  2013  (obtained  from  http://weather.uwyo.edu/ 
upperair/sounding.html).  A  full  wind  barb  represents  10  kts  and  a  half  wind  barb  represents 
5  kts  (continued). 
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Hanley  et  al.  (2016)  employed  the  United  Kingdom  Met  Office’s  Unified  Model 
with  4.4/2.2/0.5/0.2/0.1-km  horizontal  grid  spacing  nests  to  investigate  the  20 
May  2013  Moore  tornado.  They  used  no  data  assimilation  (data  assimilation  was 
applied  to  the  global  model  driving  the  outermost  nest).  Their  2.2-km  grid 
initiates  convection  at  about  the  correct  time  in  the  Oklahoma  City  area  and 
simulates  supercells,  but  convection  persists  too  long  in  this  area.  Their  200-  and 
100-m  grids  simulate  tornado-like  vortices,  but  do  so  approximately  2.5  h  later 
than  the  observed  Moore  tornado. 

Zhang  et  al.  (2015,  2016)  examined  predictability  of  the  convection  on  20  May 
2013  using  WRF-ARW.  Zhang  et  al.  (2015)  used  nested  27-,  9-,  3-,  and  1-km 
horizontal  grid  spacing  domains  to  investigate  the  impact  on  the  simulation  of 
convection  from  temporal  variations  in  initial  conditions  and  from  horizontally 
shifting  the  underlying  terrain.  They  find  that  temporal  shifting  of  initial 
conditions  generally  temporally  shifts  convection  as  expected  but  in  some  cases 
does  not  affect  convective  timing  because  the  lateral  boundary  conditions  are 
controlling  convective  initiation.  The  temporal  shifting  impacts  the  paths  of  some 
supercells  (diagnosed  via  updraft  helicity)  much  more  than  others,  suggesting  that 
the  predictability  of  some  supercells  is  larger  than  other  supercells.  Terrain 
shifting  changed  the  timing  of  convective  initiation,  as  well  as  how  the  storms 
developed.  Zhang  et  al.  (2016)  also  use  nested  27-,  9-,  3-,  and  1-km  horizontal 
grid  spacing  domains,  but  add  a  250-m  nest  to  create  perturbations  to  create  a 
60-member,  1-km  ensemble  to  investigate  predictability.  The  magnitude  of  these 
perturbations  is  unobservable  with  current  observing  platforms,  thus  the  ensemble 
represents  the  effects  of  currently  unavoidable  observational  uncertainty.  Zhang  et 
al.  (2016)  find  that  while  a  line  of  storms  is  produced  in  all  of  the  ensemble 
members,  the  details  of  individual  storms  in  that  line  differ  among  the  ensemble 
members;  this  suggests  that  the  line  of  storms  is  predictable  but  details  along  that 
line  of  storms  may  not  be  predictable  at  this  time. 

Snook  et  al.  (2016)  produced  500-m  horizontal  grid  spacing  Advanced  Regional 
Prediction  System  (ARPS)  simulations  to  investigate  numerical  weather 
prediction  forecasts  of  hail  for  this  case.  They  applied  EnKF  using  40  ensemble 
members  with  radar  reflectivity,  radar  radial  velocity,  and  surface  observations 
assimilated  every  5  min  between  1830  and  1930  UTC,  followed  by  a  90-min 
forecast  (1930-2100  UTC).  Their  results  indicated  that  for  this  case  the  model 
had  skill  at  predicting  hail. 
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4.  Methodology 


This  study  investigates  3  methods  to  produce  initial  conditions  as  well  as  the 
application  of  analysis  and  observation  nudging.  This  section  describes  the 
methodology  used  to  create  the  initial  conditions  as  well  as  the  analyses  and 
observations  that  are  nudged  toward. 

4.1  Observations  for  Use  in  Observation  Nudging 

The  observations  used  for  observation  nudging  data  assimilation  were  obtained 
from  the  Meteorological  Assimilation  Data  Ingest  System  (MADIS; 
https://madis.noaa.gov/').  The  specific  MADIS  data  sources  applied  here  are 
standard  surface,  mesonet  surface,  profiler,  radiosonde,  and  Aircraft 
Communications  Addressing  and  Reporting  System  (ACARS)  observations.  After 
simulations  were  complete,  we  discovered  that  a  software  error  resulted  in  profiler 
and  ACARS  observations  not  being  included  at  13,  19,  and  21  UTC. 

Quality  control  is  applied  using  3  resources.  First,  the  quality  control  flags 
attached  to  the  MADIS  data  set  are  used  to  remove  observations  marked  as 
problematic.  Second,  use/reject  lists  that  were  designed  for  the  Real-Time 
Mesoscale  Analysis  (De  Pondeca  et  al.  2011)  are  applied  only  to  the  surface 
mesonet  observations;  these  lists  are  designed  to  deal  with  the  data  quality  issues 
that  can  be  more  prevalent  in  the  mesonet  observations  (e.g.,  poor  siting)  than  in 
standard  surface  observations.  Third,  the  Obsgrid  program  (NCAR  2016)  is  used 
to  apply  quality  control  based  on  evaluating  observations  against  a  first-guess 
field  and  against  nearby  observations  (buddy  check).  The  first-guess  fields  we  use 
are  from  the  15  UTC  cycle  of  the  3-km  High-Resolution  Rapid  Refresh  model 
(HRRR;  Alexander  et  al.  2013;  http://rapidrefresh.noaa.gov/hrrr)  on  20  May 
2013.  Obsgrid  performs  quality  control  of  temperature,  wind,  relative  humidity, 
sea-level  pressure,  and  surface  pressure  (NCAR  2016;  Reen  2015).  In  this  study, 
the  first-guess  field  (HRRR)  is  vertically  interpolated  to  additional  pressure  levels 
to  facilitate  quality  control  of  single-level  above-surface  observations  (ACARS 
data)  and  allow  for  additional  vertical  structure  to  be  retained  in  multilevel 
observations  (e.g.,  radiosonde)  data. 

4.2  Obsgrid  Analyses 

Obsgrid  produced  an  objective  analysis  used  as  initial  conditions  for  some 
simulations.  The  analysis  combines  the  observations  in  Section  4.1  with  3-km 
HRRR  fields  from  the  15  UTC  cycle  on  20  May  2013.  Obsgrid  was  configured  here 
to  perform  multiscan  Cressman  analyses  of  temperature,  wind,  relative  humidity, 
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and  surface  pressure.  The  surface  pressure  analysis  is  used  to  determine  how  the 
surface  analyses  of  the  other  fields  are  applied  in  the  WRF  initial  conditions  based 
on  the  difference  between  the  WRF  surface  pressure  and  the  surface  pressure  of  the 
Obsgrid  analyses.  The  Cressman  analyses  for  wind  and  relative  humidity  are  flow 
dependent.  Five  scans  are  applied  using  the  default  Obsgrid  technique  of 
determining  radius  of  influence  based  on  assuming  an  average  of  325  km  between 
upper  air  observations.  This  results  in  radii  of  influence  of  520,  354,  255,  179,  and 
126  km  for  above-surface  observations.  Surface  observations  are  assumed  to  have  a 
40%  smaller  radii  of  influence,  given  the  greater  degree  of  spatial  heterogeneity  at 
the  surface;  however,  Obsgrid  limits  the  radii  of  influence  for  surface  observations 
specified  as  a  percentage  of  the  radius  of  influence  of  the  upper-air  observations  to 
be  between  4.5  and  100.0  model  grid  cells.  The  resulting  surface  radii  of  influence 
for  the  5  scans  are  100,  70, 49,  34,  and  24  km. 

4.3  vLAPS  Analyses 

For  this  case  15-min  vLAPS  analyses  were  created  using  the  3-km  HRRR  fields 
from  the  15  UTC  cycle  on  20  May  2013.  One  source  of  observations  used  in  the 
analysis  was  the  MADIS  data  set.  The  observation  types  obtained  from  the 
MADIS  data  set  included  standard  surface,  mesonet  surface,  profiler,  radiosonde, 
and  ACARS.  Pilot  observations  were  obtained  from  the  Federal  Aviation 
Administration  reporting  sky  cover,  turbulence,  wind  direction,  and  icing.  In 
addition  to  using  the  MADIS  quality  control  flag,  vLAPS  also  checks 
observations  against  the  background  field  to  determine  if  the  difference  is  large 
enough  that  the  observation  should  be  omitted  from  the  analysis.  Radial  velocity 
from  26  radars  (WSR-88D  and  Terminal  Doppler  Weather  Radar)  were  ingested 
using  the  variational  data  assimilation,  but  the  reflectivity  was  ingested  using  the 
earlier  nonvariational  LAPS.  The  LAPS  cloud  analysis  uses  GOES  imagery  data 
to  remove  clouds  where  GOES  indicates  that  clouds  are  not  present.  The  LAPS 
cloud  analysis  also  uses  radar,  surface  observations  (METARs),  and  model  first- 
guess  fields  and  uses  these  to  provide  3-D  hydrometeor  content;  the  results  of  the 
cloud  analysis  also  are  used  to  adjust  the  3-D  wind  field.  Among  the  fields 
analyzed  by  LAPS  are  temperature,  water  vapor,  winds,  and  hydrometeors. 

4.4  Creating  WRF  Input  Files 

The  process  of  creating  the  WRF-ARW  input  files  for  the  various  experiments 
described  in  Section  5  is  somewhat  involved  (Fig.  6).  In  the  following  description, 
a  letter  in  brackets  is  used  to  signify  the  component  labeled  with  that  letter  in  the 
figure.  Most  of  the  components  are  standard  WRF  Preprocessing  System  (WPS) 
software,  or  altered  versions  of  the  standard  preprocessors. 
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Fig.  6  Flowchart  illustrating  creation  of  the  WRF  input  files  for  this  study.  Colored  portions  of  the  chart  indicate  steps  required  to  create  WRF 
initial  conditions  and  boundary  conditions  of  a  certain  type  (HRRR  in  orange,  vLAPS  in  red,  and  Obsgrid  in  blue).  The  components  are  labeled  with 
letters  to  facilitate  discussion  in  the  text. 


In  the  first  step,  static  geographical  data  (e.g.,  terrain  height,  soil  type)  are 
ingested  by  Geogrid  [A]  to  create  static  geographical  data  for  the  current  domain 
(geo_em*).  The  hourly  3-km  HRRR  data  from  the  15  UTC  cycle  on  20  May  2013 
is  then  converted  from  Gridded  Binary  (GRIB)  format  into  WRF  intermediate- 
format  files  by  both  the  standard  version  of  Ungrib  [E]  and  a  “high-resolution” 
version  of  Ungrib  [B],  The  standard  version  of  Ungrib  [E]  retains  data  on  all  of 
the  pressure  levels  present  in  the  original  GRIB  file  and  does  not  add  data  on 
additional  pressure  levels.  The  high-resolution  version  of  Ungrib  [B]  differs  from 
the  standard  version  in  that  it  vertically  interpolates  the  HRRR  data  to  additional 
pressure  levels  to  allow  for  improved  quality  control  of  single-level  above-surface 
observations  and  to  better  retain  vertical  resolution  of  multilevel  observations. 
When  observations  are  quality  controlled  (at  a  later  step),  single-level  above¬ 
surface  observations  (e.g.,  aircraft  data)  are  quality  controlled  against  the  nearest 
pressure  level.  By  interpolating  the  HRRR  data  to  additional  pressure  levels,  the 
quality  control  of  these  observations  can  use  HRRR-derived  data  closer  to  the 
pressure  level  of  the  observation.  When  multilevel  observations  (e.g., 
radiosondes)  are  quality  controlled  (at  a  later  step),  they  are  vertically  interpolated 
to  the  pressure  levels  of  the  HRRR-derived  data  to  allow  quality  control  of  the 
observations  against  the  HRRR  data.  By  interpolating  the  HRRR  data  to 
additional  pressure  levels,  a  more  detailed  vertical  structure  can  be  retained  in 
multilevel  observations.  The  intermediate  files  created  by  both  versions  of  Ungrib 
([B]  and  [E])  are  ingested  by  Metgrid  ([C]  and  [F],  respectively)  to  create  HRRR 
data  gridded  to  the  WRF  domain’s  horizontal  grid  (met_em*  files). 

The  met_em*  files  originating  from  the  Metgrid  [C]  run  using  output  from  the 
high-resolution  version  of  Ungrib  [B]  are  ingested  by  Obsgrid  [D]  along  with 
MADIS  observations.  These  metem*  files  contain  the  pressure  levels  in  the 
original  HRRR  data,  but  also  include  additional  pressure  levels  onto  which  the 
HRRR  data  are  interpolated.  Obsgrid  performs  quality  control  of  the  MADIS 
observations  by  comparing  each  MADIS  observation  to  nearby  observations  and 
to  the  HRRR  data  in  the  met_em*  files  with  the  additional  pressure  levels.  The 
quality  controlled  observations  are  output  in  WRF  observation  nudging  format  as 
OBSDOMAIN*  files  ready  for  WRF  [Q]  to  ingest  for  observation  nudging. 

The  met  em*  files  originating  from  the  Metgrid  [F]  run  using  output  from  the 
non-high-resolution  version  of  Ungrib  [E]  are  used  for  creating  initial  conditions 
and  boundary  conditions.  These  met_em*  files  retain  the  pressure  levels  in  the 
original  HRRR  data.  The  subsequent  processing  of  these  met_em*  files  differs 
depending  on  which  type  of  initial  conditions  are  being  generated. 

To  create  Obsgrid-based  initial  conditions  the  met_em*  files  from  Metgrid  [F]  are 
used  as  the  first  guess  for  an  objective  analysis  created  by  Obsgrid  [H]  by 
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combining  the  HRRR-based  met  em*  files  with  the  MADIS  observations.  The 
resultant  objective  analyses  (metoa_em*)  then  have  snow-cover-related  fields 
removed  to  avoid  potential  discrepancies  with  vLAPS-based  initial  conditions  that 
do  not  have  snow  cover  as  a  field.  (Snow  cover  here  refers  to  a  yes/no  flag 
indicating  snow  cover  and  is  separate  from  snow  depth.  Snow  cover  is  not 
processed  in  the  standard  WPS/WRF  software  but  was  added  in  the  US  Army 
Research  Laboratory  [ARL]  implementation  to  allow  snow  cover  fields  to  be 
ingested  and  thus  improve  snow  characterization  in  WRF  simulations.)  The 
objective  analyses  file  for  the  initial  time  (meteoem*  from  [I])  and  the  HRRR- 
based  met_em*  with  snow  cover  removed  (met  em*  from  [G])  are  then  processed 
by  the  WRF  component  Real  to  create  initial  condition  (wrfinput*)  and  boundary 
condition  (wrfbdy*)  files  for  WRF  that  are  on  the  WRF  vertical  levels. 

The  simplest  of  the  initial  conditions  to  create  are  those  that  are  based  solely  on 
the  HRRR  output  and  thus  use  no  observations.  To  create  HRRR-based  initial 
conditions,  Real  [K]  creates  initial  condition  (wrfinput*)  and  boundary  condition 
(wrfbdy*)  files  on  the  WRF  vertical  levels  using  the  met  em*  files  from  Metgrid 
[F]  with  snow  cover  removed  [G], 

To  create  the  vLAPS-based  WRF  initial  conditions,  the  vLAPS  analyses  are 
ingested  by  Metgrid  [N]  to  create  vLAPS-based  met  em*  files.  Real  [O]  is  then 
used  to  create  initial  condition  files  (wrfinput*),  which  are  renamed 
wrfvar_output*  for  input  into  da_update_bc  [L],  The  WRF  input  and  boundary 
condition  files  created  by  Real  [K]  from  the  HRRR-based  met  em*  files  are  also 
input  to  da  update  bc  [L],  This  program  [L]  updates  the  boundary  conditions  in 
wrfbdy*  based  on  the  differences  between  the  initial  conditions  from  which  these 
boundary  conditions  were  created  (wrfinput*)  and  the  initial  conditions  that  will 
actually  be  used  (wrfvar  output*).  This  prevents  a  mismatch  near  the  boundary 
between  the  boundary  condition  and  initial  condition  file.  The  boundary  condition 
updating  program  (da_update_bc  [L])  creates  an  updated  boundary  condition  file 
(wrfbdy*)  ready  for  ingestion  by  WRF  [Q],  as  well  as  an  initial  condition  file 
wrfvar_output*  that  we  rename  wrfinput*.  This  vLAPS-based  initial  condition 
file  does  not  contain  the  HRRR  soil  moisture  and  soil  temperature  fields,  nor  the 
netCDF  file  “global  attributes”.  Therefore,  we  use  the  tool  NetCDF  Operators 
(http://nco.sourceforge.net/)  component  neks  [M]  to  take  the  soil  moisture,  soil 
temperature,  and  global  attributes  from  the  HRRR-based  initial  conditions  (the 
wrfinput*  created  by  Real  [K])  and  place  them  in  the  vLAPS-based  initial 
conditions  (wrfinput*,  which  is  the  wrfvar_output*  files  output  by  da_update_bc 
[L]).  The  resultant  initial  condition  file  is  the  vLAPS  wrfinput*  file  used  by  WRF 

[Q]- 


Approved  for  public  release;  distribution  is  unlimited. 


28 


In  addition  to  the  observation  nudging  file  (OBS  DOMAIN*)  and  the  initial  and 
boundary  conditions  (wrfinput*  and  wrfbdy*),  we  must  also  produce  surface  and 
3-D  analysis  nudging  files.  These  are  created  from  the  vLAPS-based  metem* 
files  created  by  Metgrid  [N],  Applying  Real  [O]  to  these  files  creates  the  3-D 
analysis  nudging  files  (wrffdda*).  Using  Obsgrid  [P]  on  these  met_em*  files 
without  providing  observations  to  Obsgrid  [P]  creates  vLAPS-based  surface 
analysis  nudging  files  (wrfsfdda*)  ready  for  ingestion  by  WRF  [Q], 


5.  Experimental  Design 


A  series  of  experiments  (Table  1;  Fig.  7)  were  carried  out  to  investigate  the 
potential  value  of  combining  the  vLAPS  and  nudging  data  assimilation 
techniques.  The  experiments  differed  in  their  initial  conditions  (HRRR,  vLAPS, 
or  Obsgrid),  the  length  of  pre-forecast  applied  (0,  1,  or  3  h),  and  whether 
observation  or  analysis  nudging  was  applied.  The  pre-forecast  is  used  to  refer  to  a 
period  at  the  beginning  of  the  model  integration  during  which  1)  observations  are 
assimilated  or  2)  the  model  is  assumed  to  be  spinning  up.  Assimilation  of 
observations  taken  during  the  pre-forecast  may  extend  into  the  beginning  of  the 
free  forecast,  but  no  observations  valid  during  the  free  forecast  should  be 
assimilated.  However,  the  analyses  being  analysis  nudged  toward  here  include 
observations  valid  during  the  beginning  of  the  forecast  period  (first  30  min  for 
Obsgrid  analyses  and  15  min  for  vLAPS  analyses).  The  experiments’  names  start 
with  text  indicating  the  initial  condition  source  (“HRRR”,  “VLAPS”,  or 
“OBSGRID”),  then  the  length  of  the  pre-forecast  in  hours  (single  integer), 
followed  by  “O”  if  observation  nudging  is  applied,  and  “A”  if  analysis  nudging  is 
applied. 


Table  1  Experimental  design 


Assimilation 

Name 

Initial 

Condition 

Source 

Pre¬ 
forecast 
Length  (h) 

Nudging 

Analysis  Obs 

None 

HRRRO 

HRRR 

0 

N 

N 

None 

HRRR3 

HRRR 

3 

N 

N 

vLAPS 

VLAPSO 

VLAPS 

0 

N 

N 

vLAPS 

VLAPS  1 

VLAPS 

1 

N 

N 

vLAPS 

VLAPS3 

VLAPS 

3 

N 

N 

Ob  sgrid/Nudging 

0BSGRID30 

OBSGRID 

3 

N 

Y 

Partial  hybrid 

VLAPS30 

VLAPS 

3 

N 

Y 

Partial  hybrid 

VLAPS3A 

VLAPS 

3 

Y 

N 

Full  hybrid 

VLAPS3AO 

VLAPS 

3 

Y 

Y 
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Fig.  7  Schematic  showing  experimental  design,  to  indicates  the  end  of  the  pre-forecast  and 
the  beginning  of  the  forecast,  and  X  will  vary  among  model  cycles,  but  all  model  cycles  end 
at  00  UTC. 


The  first  experiment  (HRRRO)  serves  as  a  control  experiment.  It  uses  HRRR  to 
determine  initial  and  boundary  conditions  and  uses  neither  a  pre-forecast  nor 
nudging.  HRRR3  is  similar  but  includes  a  3-h  model  spin-up  time  period  to 
investigate  the  potential  benefit  of  a  3-h  pre-forecast  even  if  no  data  assimilation 
is  applied.  This  will  also  facilitate  evaluation  of  other  experiments  and  whether 
differences  between  experiments  are  due  to  a  3-h  data  assimilation  period  or 
simply  due  to  starting  the  model  integration  3  h  earlier. 
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The  second  set  of  experiments  (VLAPSO,  VLAPS1,  and  VLAPS3)  uses  the 
vLAPS  analysis  as  the  initial  condition  and  uses  either  no  pre-forecast  spin-up 
period  (VLAPSO),  a  1-h  pre-forecast  spin-up  period  (VLAPS  1),  or  a  3-h  pre¬ 
forecast  spin-up  period  (VLAPS3).  The  VLAPSO  experiment  most  closely 
matches  applications  of  the  vLAPS  technique  in  past  studies,  but  the  other  2 
experiments  are  added  here  to  determine  the  potential  benefit  and  preferred  length 
of  a  spin-up  period;  VLAPS3  also  facilitates  comparison  with  nudging 
experiments  wherein  nudging  is  applied  over  3  h. 

0BSGRID30  uses  an  objective  analysis  created  by  Obsgrid  (see  Section  4.2)  as 
initial  conditions  for  WRF  and  then  applies  observation  nudging  also  created  via 
Obsgrid  (see  Section  4.1)  over  the  first  3  h.  0BSGRID30  matches  the  data 
assimilation  methodology  commonly  used  at  ARL. 

The  final  3  experiments  combine  vLAPS  analysis  with  nudging.  VLAPS30  uses 
vLAPS  to  specify  the  initial  conditions  and  then  applies  observation  nudging  over 
the  first  3  h.  VLAPS3A  also  uses  vLAPS  to  specify  the  initial  conditions,  but  also 
analysis  nudges  toward  15-min  vLAPS  analyses  during  the  first  3  h  (it  does  not 
apply  any  observation  nudging).  VLAPS3AO  is  the  full  hybrid  experiment 
wherein  vLAPS  analysis  are  used  to  specify  the  initial  conditions  and  are  nudged 
toward  during  the  first  3  h,  but  the  model  solution  is  also  observation  nudged 
toward  observations  during  these  same  first  3  h. 

Multiple  cycles  were  completed  for  each  experiment.  The  0-h  forecast  (i.e.,  the 
beginning  of  the  forecast  period)  for  the  first  cycle  in  each  experiment  is  18  UTC. 
However,  the  model  start  time  of  the  first  cycle  is  either  15,  17,  or  18  UTC 
depending  on  the  pre-forecast  length  for  that  experiment  (Table  1).  Each  cycle  of 
each  experiment  ends  at  00  UTC,  which  means  that  each  experiment  has  a  6-h 
forecast  for  the  first  cycle. 

The  beginning  of  the  forecast  for  the  next  cycle  of  each  experiment  is  at  19  UTC, 
and  the  ending  of  the  forecast  is  at  00  UTC;  thus,  each  experiment  has  a  5-h 
forecast.  Each  cycle  used  only  HRRR  output  from  HRRR  integration  that  begins 
at  15  UTC.  The  cycles  are  independent,  thus  no  data  from  the  WRF  integration  in 
the  first  cycle  are  used  in  the  WRF  integration  in  the  second  cycle.  The  cycles 
continue  until  the  23  UTC  cycle,  which  produces  a  1-h  forecast. 

Each  cycle  is  referred  to  by  the  beginning  of  the  forecast  period  of  the  cycle  (to  in 
Fig.  7).  For  example,  the  VLAPS3  cycle  that  starts  at  16  UTC  is  referred  to  as 
VLAPS3i9  since  the  first  3  h  are  pre-forecast  (Table  1)  and  thus  the  forecast  starts 
at  19  UTC. 
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6.  Preliminary  Results 


Currently,  only  a  preliminary  analysis  of  the  results  of  these  experiments  has  been 
completed.  This  technical  report  is  intended  primarily  to  describe  the  methods 
employed  in  this  study  and  provide  some  preliminary  results,  with  a  subsequent 
technical  report  providing  a  more  detailed  analysis  of  the  results. 

6.1  Experiment  Comparison  at  Time  of  Newcastle-Moore 
Tornado 

At  2015  UTC,  the  Newcastle-Moore  tornado  was  at  or  near  its  peak  strength.  As 
part  of  the  preliminary  analysis  of  the  experiment  results,  we  investigated 
differences  between  the  experiments  at  this  time.  The  observed  radar  composite 
base  reflectivity  at  2015  UTC  (Fig.  8)  shows  a  broken  southwest-northeast 
oriented  line  of  strong  convection  passing  through  central  Oklahoma  and 
southeastern  Kansas.  The  Newcastle-Moore  tornado  (black  “X”  in  the  white 
square)  is  near  the  rear  of  this  line.  A  weak,  similarly  oriented  line  of  precipitation 
is  located  over  southcentral  Kansas.  There  is  also  scattered  precipitation  in 
western  Kansas  and  Oklahoma,  and  in  the  portion  of  southeastern  Colorado  in  the 
domain. 
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Fig.  8  Composite  base  reflectivity  at  2015  UTC.  The  black  X  inside  the  white  square 
(near  the  end  of  the  white  arrow)  indicates  the  location  of  the  Newcastle-Moore  tornado  at 
2015  UTC.  Radar  reflectivity  obtained  from  Iowa  Environmental  Mesonet 
(http://mesonet.agron.iastate.edu/GIS/apps/rview/warnings.phtml)  and  cropped  to  the 
approximate  extent  of  the  WRF-ARW  model  domain. 

6.1.1  Comparison  of  Hybrid  to  Component  Techniques 

The  first  4  experiments  examined  represent  a  control  experiment,  the  2  components 
of  the  hybrid,  and  the  full  hybrid  experiment.  Namely,  HRRR3  is  chosen  as  the 
control  since  it  contains  no  data  assimilation  but  has  a  3-h  pre-forecast  to  allow 
convection  to  spin  up.  VLAPSO  represents  the  vLAPS  component  of  the  hybrid, 
since  it  represents  how  vLAPS  is  usually  applied  (i.e.,  with  no  pre-forecast  spin- 
up).  0BSGRID30  represents  the  nudging  component,  with  the  initial  conditions 
enhanced  using  an  objective  analysis  created  by  Obsgrid  and  3  h  of  observation 
nudging  applied.  VLAPS3AO  is  the  full  hybrid  experiment,  incorporating  vLAPS 
as  the  initial  conditions,  analysis  nudging  to  15 -min  vLAPS  analyses,  and 
observation  nudging.  We  first  examine  the  results  from  the  experiments  whose  0-h 

Approved  for  public  release;  distribution  is  unlimited. 

33 


forecast  is  at  18  UTC  (i.e.,  HRRR3is,  VLAPSOis,  0BSGRID30is,  and 
VLAPS3AOi8). 


The  model-simulated  lowest  model  level  radar  reflectivity  of  all  4  experiments  at 
2015  UTC  (Fig.  9)  generally  reproduces  the  overall  structures  in  the  observed 
composite  base  reflectivity  (Fig.  8).  In  particular,  all  4  experiments  show  the  main 
line  of  convection  along  the  proper  alignment,  scattered  precipitation  in  the 
northwestern  comer  of  the  domain,  and  some  precipitation  in  the  northcentral 
portion  of  the  domain.  However,  there  are  significant  differences  in  the  details  of 
predicted  precipitation. 
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Fig.  9  WRF-ARW  lowest  model  level  radar  reflectivity  at  2015  UTC  on  20  May  2013  for 
a)  HRRR3i8,  b)  VLAPSOis,  c)  0BSGRID30is,  and  d)  VLAPS3AOis.  The  location  of  the 
Newcastle-Moore  tornado  at  2015  UTC  is  indicated  by  the  black  “X”  within  the  white 
square  (in  the  same  location  as  marked  by  the  arrow  in  Fig.  8). 
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There  is  considerable  variation  among  the  experiments  in  the  prediction  of  the 
area  of  light  precipitation  observed  over  southcentral  Kansas  (Fig.  8;  top  center  of 
plot).  HRRR3i8  (Fig.  9a)  and  even  more  so  0BSGRID30is  (Fig.  9c)  underpredict 
the  coverage  of  this  feature.  VLAPSOis  (Fig.  9b)  predicts  2  roughly  parallel  lines 
of  convection,  with  the  leading  line  much  stronger  than  the  observed  structure, 
and  the  trailing  line  more  scattered  than  observed.  While  VLAPS3AOi8  (Fig.  9d) 
predicts  higher  reflectivity  than  observed  and  with  less  coverage  than  observed,  it 
appears  to  best  predict  the  coverage  of  this  precipitation. 

The  portion  of  the  main  convective  line  over  southeastern  Kansas  (Fig.  8)  shows 
substantial  differences  between  the  experiments.  HRRR3i8  (Fig.  9a)  and 
VLAPSOis  (Fig.  9b)  predict  convection  that  is  much  more  scattered  than  the 
continuous  line  observed,  although  VLAPSOis  appears  to  perform  better  than 
HRRR3i8  as  it  indicates  greater  coverage  than  HRRR3i8.  0BSGRID30is 
(Fig.  9c)  predicts  stronger  convection  with  much  more  coverage  than  HRRRis  and 
VLAPSO  is  but  does  not  simulate  the  linear  nature  of  the  convection  as  well  as 
VLAPS3AOi8  (Fig.  9d). 

Over  Oklahoma,  the  predictions  of  the  main  convective  line  also  differ  among  the 
experiments.  Convection  extends  too  far  into  the  northeastern  comer  of  Oklahoma 
in  0BSGRID30i8  (compare  Fig.  9c  to  Fig.  8)  and  VLAPS3AOis  (Fig.  9d). 
HRRR3i8  (Fig.  9a)  and  VLAPSOis  (Fig.  9b)  do  not  seem  to  have  this  issue.  Also, 
VLAPS3AO  is  also  appears  to  better  represent  the  discrete  nature  of  the 
convection  in  the  main  convective  line  than  the  other  experiment. 

All  4  experiments  indicate  moist  convective  elements  on  the  western  edge  of 
Oklahoma,  but  radar  indicates  that  this  area  was  free  of  precipitation  at  this  time 
(Fig.  8).  In  this  case,  0BSGRID30  is  (Fig.  9c)  predicts  weaker  convection  and 
with  less  coverage  than  any  of  the  other  experiments,  while  VLAPS3AOis 
(Fig.  9d)  appears  to  predict  less  coverage  of  stronger  convection  (>30  dBZ)  than 
the  other  2  experiments  (HRRR3is  and  VLAPSOis). 

While  all  4  experiments  predict  a  relatively  small  area  of  convection  in  the  far 
south  of  the  central  portion  of  the  domain  (Fig.  9)  consistent  with  observations 
(Fig.  8),  some  predict  convection  a  couple  counties  north  of  this  location,  which  is 
inconsistent  with  observed  radar.  VLAPS3AOis  (Fig.  9d)  predicts  this  convection 
to  be  the  most  widespread  while  0BSGRID30is  (Fig.  9c)  also  predicts  this  area 
of  convection  but  forecast  that  it  would  not  be  as  strong  as  indicated  by 
VLAPS3AOi8. 

These  4  experiments  vary  in  the  predictions  near  the  observed  location  of  the 
Newcastle-Moore  tornado  at  2015  UTC.  0BSGRID30i8  (Fig.  9c)  shows 
convection  just  northwest  of  this  location  while  HRRR3is  (Fig.  9a)  and 
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VLAPS3AOi8  (Fig.  9d)  show  convection  north  of  this  location.  VLAPSOis 
(Fig.  9b)  indicate  convection  north  of  this  location.  VLAPSOis  (Fig.  9b)  is  the 
only  experiment  that  predicts  convection  at  the  location  of  the  Newcastle-Moore 
tornado  at  this  time. 

6.1.2  Comparison  between  Varieties  of  Component  Techniques 

While  the  previous  section  explored  the  difference  between  the  full  hybrid  versus 
its  components,  here  we  investigate  differences  between  varieties  of  the 
component  techniques  that  differ  in  the  length  of  the  pre-forecast  or  on  what  types 
of  nudging  are  applied.  In  this  section,  we  again  solely  examine  the  cycle  with  an 
18  UTC  0-h  forecast.  First,  we  compare  HRRROis  versus  HRRR3is,  then 
VLAPSOis  versus  VLAPSI18  and  VLAPS3is,  and  finally  VLAPS3AOi8  versus 
VLAPS3Ai8  and  VLAPS30i8. 

Comparing  the  2  experiments  that  use  neither  vLAPS  analyses  nor  nudging 
(HRRRO  is  and  HRRR3i8)  reveals  the  potential  benefit  of  a  3-h  period  for  the 
model  spin-up.  The  experiments  only  differ  in  that  HRRROis  starts  at  18  UTC, 
whereas  HRRR3is  starts  at  15  UTC  and  thus  contains  3  h  (15-18  UTC)  during 
which  model  spin-up  can  occur.  The  lowest  model  level  model-simulated 
reflectivity  in  HRRROis  (Fig.  10b)  is  very  similar  to  HRRR3is  (Fig.  10a)  at  2015 
UTC.  Without  the  3-h  spin-up,  HRRROis  still  reproduces  the  same  overall 
structure  of  precipitation  as  HRRR3is.  However,  for  the  time  under  evaluation 
(2015  UTC),  even  though  HRRROis  does  not  have  the  15-18  UTC  time  for  spin- 
up,  it  has  had  just  over  2  h  to  spin  up  the  convection.  HRRROis  increases  the 
coverage  of  the  precipitation  halfway  across  the  northern  edge  of  the  domain  that 
HRRR3  is  underpredicted  but  shows  the  precipitation  as  more  intense  than  radar 
indicated  (Fig.  10c)  and  does  not  extend  the  convection  nearly  far  enough  to  the 
southwest.  Contrary  to  observations  and  HRRR3is,  HRRROis  produces  scattered 
convection  between  this  element  (the  precipitation  halfway  across  the  northern 
edge  of  the  domain)  and  the  widespread  scattered  convection  in  the  northwest 
comer  of  the  domain.  Both  experiments  fail  to  reproduce  the  line  of  convection 
across  southeastern  Kansas.  In  the  area  of  the  Newcastle-Moore  tornado,  moist 
convection  is  somewhat  farther  away  to  the  north  but  much  closer  to  the  south  in 
HRRROis  as  compared  to  HRRR3is.  Overall,  it  is  not  clear  whether  HRRROis  or 
HRRR3  is  performs  better  at  this  time.  The  similarities  between  HRRROis  and 
HRRR3i8  are  not  entirely  unexpected  since  both  are  initialized  using  output  from 
the  same  3 -km  HRRR  simulation  and  since  both  experiments  are  using  physics 
options  very  similar  to  those  used  in  the  HRRR  simulation. 
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Fig.  10  WRF-ARW  lowest  model  level  radar  reflectivity  at  2015  UTC  on  20  May  2013  for 
a)  HRRR3i8,  b)  HRRROis,  and  c)  observed  base  reflectivity  composite.  The  location  of  the 
Newcastle-Moore  tornado  at  2015  UTC  is  indicated  by  the  black  “X”  within  the  white 
square. 

The  experiments  VLAPSOis,  VLAPSlis,  and  VLAPS3is  assist  in  evaluating 
whether  using  a  more  recent  vLAPS  analysis  for  initial  conditions  overwhelms 
whatever  benefits  are  accrued  via  additional  model  spin-up  time.  At  the  2015 
UTC  time  shown,  experiments  VLAPSOis,  VLAPSlis,  and  VLAPS3is  have 
integrated  2.25,  3.25,  and  5.25  h,  respectively,  after  accounting  for  the  varying 
pre-forecast  times  assigned  (0,  1,  and  3  h)  and  the  2.25  h  between  the  1800  UTC 
start  of  this  cycle  and  the  2015  UTC  time  being  discussed.  The  weak  line  of 
precipitation  in  central  Kansas  (extending  toward  the  southwest  from  the  north- 
central  edge  of  the  observed  radar;  Fig.  lid)  was  overpredicted  in  VLAPSOis 
(Fig.  11a),  while  adding  an  hour  pre-forecast  (VLAPSlis;  Fig.  lib)  significantly 
weakened  this  precipitation,  an  additional  2  h  of  pre-forecast  (VLAPS3is; 
Fig.  lie)  results  in  a  structure  that  better  resembles  that  observed.  The  line  of 
convection  in  southeastern  Kansas  becomes  better  organized  with  increased  pre- 
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forecast  time  and  more  closely  matches  the  line  observed.  The  southern  end  of  the 
main  line  of  convection  varies  among  the  experiments,  but  it  is  difficult  to 
determine  which  performs  best.  The  spatial  extent  of  the  spurious  convection  in 
the  northwestern  portion  of  Oklahoma  appears  to  be  minimized  with  a  3-h  pre¬ 
forecast,  although  a  0-h  pre-forecast  outperforms  a  1-h  pre-forecast.  Overall,  it 
appears  that  the  3-h  pre-forecast  experiment  performs  best  in  terms  of  simulated 
reflectivity  at  2015  UTC  for  the  18  UTC  cycle.  This  suggests  the  potential 
importance  of  allowing  at  least  some  additional  model  spin-up  time  following  a 
vLAPS  “hot  start”  analysis. 
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Fig.  11  WRF-ARW  lowest  model  level  radar  reflectivity  at  2015  UTC  on  20  May  2013  for 
a)  VLAPSOis,  b)  VLAPSlis,  c)  VLAPS3i8,  and  d)  observed  base  reflectivity  composite.  The 
location  of  the  Newcastle-Moore  tornado  at  2015  UTC  is  indicated  by  the  black  “X”  within 
the  white  square. 

The  next  comparison  evaluates  the  relative  contributions  of  the  analysis  nudging 
versus  the  observation  nudging  in  the  hybrid  data  assimilation  using  3 
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experiments  that  all  are  initialized  with  the  vLAPS  analysis  and  include  a  3-h  pre¬ 
forecast  time  period  with  data  assimilation.  In  addition  to  being  initialized  with 
the  vLAPS  analysis,  VLAPS3Ais  also  analysis  nudges  toward  15-min  vLAPS 
analyses  during  the  3-h  pre-forecast.  VLAPS30i8  does  not  perform  analysis 
nudging  but  includes  3  h  of  observation  nudging,  while  the  full  hybrid  experiment 
VLAPS3AOi8  performs  both.  The  full  hybrid  experiment  (VLAPS3AOi8; 
Fig.  12a)  seems  to  best  reproduce  the  northern  portion  of  the  main  line  of 
convection  (north  of  Oklahoma),  while  it  and  the  other  analysis-nudging 
experiment  (VLAPS3Ais;  Fig.  12b)  both  seem  to  best  reproduce  the  linear 
structure  across  central  Kansas  (top-center  of  plots)  compared  to  the  nonanalysis- 
nudging  experiment  (VLAPS30is;  Fig.  12c).  However,  the  nonanalysis  nudging 
experiment  (VLAPS30i8;  Fig.  12c)  best  minimizes  the  erroneously  predicted 
precipitation  in  the  western  edge  of  Oklahoma.  Just  south  of  the  southern  end  of 
the  main  convective  line,  VLAPS3AOis  produces  a  convective  element  not 
observed  (Fig.  12d).  While  VLAPS3Ais  (Fig.  12b)  and  VLAPS30is  (Fig.  12c)  do 
not  predict  an  area  of  convection  as  large  as  that  predicted  by  VLAPS3AOis,  they 
both  appear  to  produce  smaller  convective  cells  that  appear  to  be  erroneous. 
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Fig.  12  WRF-ARW  lowest  model  level  radar  reflectivity  at  2015  UTC  on  20  May  2013  for 
a)  VLAPS3AOi8,  b)  VLAPS3Ais,  c)  VLAPS30is,  and  d)  observed  base  reflectivity 
composite.  The  location  of  the  Newcastle-Moore  tornado  at  2015  UTC  is  indicated  by  the 
black  “X”  within  the  white  square. 

In  the  region  near  Moore,  VLAPS3AOi8  (Fig.  12a)  and  VLAPS3Ai8  (Fig.  12b) 
indicate  strong  convection  to  the  southeast  (stronger  in  VLAPS3AOis),  but  only 
weaker  returns  in  the  immediate  vicinity  of  Moore  (also  stronger  in 
VLAPS3AOi8).  VLAPS30i8  differs  in  that  it  indicates  convection  either  over  or 
very  close  to  Moore.  Although  there  are  differences  between  the  vLAPS 
experiments  with  nudging,  it  is  not  clear  which  experiment  performs  best  overall 
at  2015  UTC. 

6.1.3  Comparison  of  Different  Cycles 

Generally,  one  might  expect  that  forecasts  initialized  closer  in  time  to  an  event 
would  perform  better  than  forecasts  initialized  farther  in  the  past,  with  the  caveat 
that  sufficient  time  must  be  available  between  the  model  start  and  the  event  for  the 
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model  to  spin-up.  This  is  because  the  longer  the  model  integrates  without  the 
benefit  of  observations  to  constrain  the  solution  to  reality,  the  longer  errors 
introduced  by  the  initial  conditions  have  to  grow,  and  the  longer  conditions 
influenced  by  the  lateral  boundary  conditions  have  to  advect  to  the  area  of  interest. 
Therefore,  more  precisely,  assuming  sufficient  model  spin-up  time  has  already 
occurred,  one  would  expect  that  over  short-term  forecasts  the  longer  the  model 
integrates  without  being  constrained  by  observations,  the  larger  the  errors  in  the 
solution. 

The  2015  UTC  model-simulated  radar  reflectivity  from  3  HRRR3  cycles  are 
compared  to  the  observed  radar  reflectivity  in  Fig.  13.  Figure  13a  shows 
HRRR3i8  after  3-h  of  model  spin-up  (15-18  UTC)  and  2.25  h  of  subsequent 
model  integration,  Fig.  13b  shows  HRRR3i9  after  3  h  of  model  spin-up  (16-19 
UTC)  and  1.25  h  of  subsequent  model  integration,  and  Fig.  12c  shows  HRRR320 
after  3  h  of  model  spin-up  (17-20  UTC)  and  0.25  h  of  subsequent  model 
integration.  For  HRRR3  these  3  cycles  generally  show  fairly  similar  solutions  at 
2015  UTC.  All  3  simulations  show  the  main  line  of  observed  convection 
(Fig.  13d),  although  much  too  weak  in  Kansas,  scattered  convection  in  the 
northwestern  comer  of  the  domain,  a  too-weak  short  line  in  the  northcentral 
portion  of  the  domain,  and  all  inaccurately  place  convection  in  western 
Oklahoma.  None  of  the  cycles  place  strong  convection  over  Moore  at  this  time.  It 
is  important  to  note  that  HRRR3  does  not  do  any  assimilation  and  that  all  of  the 
cycles  of  HRRR3  are  driven  by  output  from  the  same  HRRR  model  forecast. 
Therefore,  the  HRRR3  cycle  that  integrates  the  shortest  before  the  2015  UTC 
time  examined  here  (HRRR320)  does  not  include  any  newer  observations  than  the 
HRRR3  cycle  that  integrates  longest  before  2015  UTC  (HRRR3is).  Thus,  it  is  not 
surprising  that  more  recent  cycles  do  not  appear  to  improve  the  forecast  at 
2015  UTC.  HRRR3  does  serve  an  important  purpose  in  the  comparison  here  in 
that  it  provides  context  for  evaluating  other  experiments  since  HRRR3  indicates 
that  change  in  the  solution  at  2015  UTC  is  caused  solely  by  initializing  using 
HRRR  at  different  times. 
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Fig.  13  WRF-ARW  lowest  model  level  radar  reflectivity  at  2015  UTC  on  20  May  2013  for 
a)  HRRR3i8,  b)  HRRR3i9,  c)  HRRR320,  and  d)  observed  base  reflectivity  composite.  The 
location  of  the  Newcastle-Moore  tornado  at  2015  UTC  is  indicated  by  the  black  “X”  within 
the  white  square. 

Unlike  the  HRRR3  simulations,  in  VLAPS3  later  cycles  include  newer 
observations  and  thus  hopefully  provide  better  forecasts.  The  2015  UTC  forecast 
in  Fig.  14a  for  VLAPS3is  is  based  on  a  WRF  simulation  started  at  15  UTC  with  a 
vLAPS  analysis,  integrated  3  h  (15-18  UTC)  to  allow  the  model  solution  to  spin 
up,  and  then  integrated  another  2.25  h  to  2015  UTC.  VLAPS3i9  (Fig.  14b)  starts 
1  h  later,  thus  the  initial  conditions  include  observations  that  are  1  h  more  recent 
than  VLAPS3i8,  and  similarly  VLAPS320  (Fig.  14c)  includes  observations  2  h 
more  recent  than  VLAPS3i8. 
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Fig.  14  WRF-ARW  lowest  model  level  radar  reflectivity  at  2015  UTC  on  20  May  2013  for 
a)  VLAPS3i8,  b)  VLAPS3i9,  c)  VLAPS320,  and  d)  observed  base  reflectivity  composite.  The 
location  of  the  Newcastle-Moore  tornado  at  2015  UTC  is  indicated  by  the  black  “X”  within 
the  white  square. 

Clear  variations  are  seen  among  the  3  cycles  of  YLAPS3  discussed  here.  The 
location  of  the  weak  line  of  precipitation  in  the  northcentral  portion  of  the  domain 
(central  Kansas)  is  predicted  well  in  VLAPS3i8  (Fig.  14a),  while  the  coverage  is 
lacking  and  the  intensity  is  larger  than  observed  (Fig.  14d).  However,  starting 
with  1-h  newer  vLAPS  analyses  (16  UTC  analysis  in  VLAPS3i9;  Fig.  14b)  nearly 
removes  this  line  and  results  in  scattered  precipitation  just  west  of  this,  both 
contrary  to  observations.  Moving  forward  to  using  the  17  UTC  vLAPS  analysis 
(VLAPS32o;  Fig.  14c)  allows  this  line  to  be  re-established  but  results  in  a  series  of 
parallel  lines  of  precipitation.  In  terms  of  this  line  of  precipitation  initializing 
WRF  with  an  earlier  vLAPS  analysis  (15  UTC)  and  allowing  WRF  a  longer  spin- 
up  time  results  in  the  forecast  that  best  matches  the  observed  nature  of  this 
precipitation.  The  northeastern  portion  of  the  main  line  of  convection  least 
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resembles  the  continuous  line  of  precipitation  observed  in  the  experiment  using 
the  newest  vLAPS  analysis  (VLAPS320;  Fig.  14c).  The  erroneous  area  of 
precipitation  in  westernmost  Oklahoma  appears  to  be  strongest  in  the  cycle  using 
the  newest  vLAPS  analysis  (VLAPS320;  Fig.  14c).  In  terms  of  the  location  of  the 
Newcastle-Moore  tornado,  all  3  experiments  indicate  this  location  is  on  the 
backside  of  the  line  of  convection,  which  is  also  indicated  by  radar.  VLAPS3is 
indicates  precipitation  at  this  location,  while  VLAPS3i9  and  VLAPS320  show  a 
line  of  convection  stretching  southwest  from  the  main  line  of  convection  to  just 
north  of  this  location  (closer  to  Newcastle-Moore  in  VLAPS320  than  in 
VLAPS3i9).  Radar  (Fig.  14d)  indicates  a  convective  element  aligned  generally 
west-east,  with  the  Newcastle-Moore  tornado  near  the  western  edge  of  this 
element.  Thus,  there  are  some  similarities  between  the  VLAPS3  and  radar- 
observed  representations  of  the  convection  in  the  region  of  the  Newcastle-Moore 
tornado. 

Both  VLAPSO  and  YLAPS3  are  initialized  with  vLAPS  analyses,  but  VLAPSO 
assigns  the  0-h  forecast  time  to  the  model  start  time,  while  VLAPS3  assigns  the 
0-h  forecast  time  3  h  after  the  model  start  time  to  allow  time  for  model  spin-up. 
Thus,  the  VLAPS320  forecast  at  2015  UTC  (Fig.  14c)  is  initialized  with  a  vLAPS 
analysis  at  17  UTC  and  integrated  3.25  h,  while  the  VLAPSO  is  forecast  at 
2015  UTC  (Fig.  15a)  is  initialized  with  a  vLAPS  analysis  at  18  UTC  and 
integrated  2.25  h.  Therefore,  Figs.  14  and  15  provide  the  2015  UTC  forecast  by 
WRF  simulations  initialized  with  vLAPS  analysis  at  15  (Fig.  14a),  16  (Fig.  14b), 
17  (Fig.  14c),  18  (Fig.  15a),  19  (Fig.  15b),  and  20  UTC  (Fig.  15c). 

Encouragingly,  the  VLAPSO20  2015  UTC  forecast  (Fig.  15c)  relatively  closely 
matches  the  radar’s  observation  at  this  time  (Fig.  15d).  Since  the  vLAPS  analyses 
used  to  initialize  the  VLAPSO  experiments  incorporate  radar  observations,  one 
would  expect  the  very  short-term  model  forecast  to  closely  resemble  the  radar 
observations  if  the  vLAPS  analysis  is  sufficiently  in  balance  with  the  model’s 
dynamics  and  physics.  If  the  vLAPS  analysis  were  substantially  out  of  balance 
with  the  WRF  model’s  dynamics  and  physics,  the  WRF  model  solution  might 
quickly  diverge  from  observations  as  large  adjustments  occur  to  bring  the  solution 
in  balance  with  the  model’s  equations.  Thus,  the  similarity  here  between  model 
and  radar  suggests  that  the  vLAPS  analyses  may  be  in  fairly  good  balance  with 
WRF  physics  and  dynamics.  However,  significant  further  analysis  would  be 
necessary  to  fully  address  this  issue. 
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Fig.  15  WRF-ARW  lowest  model  level  radar  reflectivity  at  2015  UTC  on  20  May  2013  for 
a)  VLAPSOis,  b)  VLAPSO19,  c)  VLAPSO20,  and  d)  observed  base  reflectivity  composite.  The 
location  of  the  Newcastle-Moore  tornado  at  2015  UTC  is  indicated  by  the  black  “X”  within 
the  white  square. 

As  the  analysis  time  approaches  the  forecast  time,  the  secondary  line  of 
convection,  the  one  in  the  northcentral  portion  of  the  domain,  is  better  predicted. 
However,  it  appears  to  be  overpredicted  even  in  VLAPSO20  (Fig.  15c).  The 
northeast  portion  of  the  main  line  of  convection  is  not  well  simulated  in  either 
VLAPSOis  (Fig.  15a)  or  VLAPSO19  (Fig.  15b),  but  VLAPSO20  does  much  better. 
The  west-to-east  oriented  convection  that  spawned  the  Newcastle-Moore  tornado 
is  not  well-simulated  in  any  of  the  cycles,  even  in  VLAPSO20.  The  erroneous  area 
of  convection  in  far  western  Oklahoma  is  predicted  in  both  VLAPSOis  and 
VLAPS19,  but  is  omitted  in  VLAPS20;  however,  it  may  be  that  VLAPS20  has  not 
had  sufficient  time  to  spin  up  this  convection  and  thus  will  predict  this  convection 
delayed  from  the  other  simulations.  Contrary  to  the  improvement  seen  with 
decreasing  lead  time  in  most  convective  elements,  the  southernmost  observed 
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convective  element  is  predicted  by  VLAPSOis,  but  not  by  VLAPSO19  or 
VLAPSO20.  Overall,  the  experiment  with  the  most  recent  vLAPS  analysis  appears 
to  perform  best  (VLAPS020),  but  since  the  model  has  only  had  15  min  of 
integration,  it  is  unclear  whether  this  benefit  persists  past  the  beginning  of  the 
model  forecast. 

Because  of  differences  in  the  data  assimilation  techniques  used  by  0BSGRID30 
and  VLAPS3/VLAPS0,  care  must  be  taken  in  interpreting  comparisons  between 
0BSGRID30  and  VLAPS3/VLAPS0.  0BSGRID30  has  the  advantage  over 
YLAPS3  in  that  in  addition  to  incorporating  observations  into  its  initial  conditions 
(via  the  objective  analysis  produced  by  Obsgrid)  it  also  assimilates  observations 
through  the  first  3  h  of  the  model  integration.  Thus,  the  newest  observations 
assimilated  into  0BSGRID30i8  are  at  18  UTC  while  the  newest  observations 
assimilated  into  VLsAPS3is  are  at  15  UTC.  While  VLAPSO18  is  more  comparable 
to  0BSGRID30i8  in  that  both  may  be  using  18  UTC  observations, 
0BSGRID30  is  has  the  advantage  of  3  h  of  model  spin-up  and  application  of 
observations  to  the  model  solution  over  3  h  rather  than  at  only  the  initial  time. 
Because  of  these  reasons,  one  must  be  careful  not  to  conclude  that  Obsgrid 
analyses  are  better  than  vLAPS  analysis  solely  due  to  comparisons  between 
Obsgrid-initialized  and  vLAPS-initialized  experiments  that  started  at  the  same 
time  (e.g.,  0BSGRID30i8  and  VLAPS3is).  Ultimately,  it  is  the  combination  of 
vLAPS  analyses  and  the  nudging  technique  that  we  hypothesize  will  create  the 
best  forecast. 

The  0BSGRID30  solution  varies  among  the  model  start  times  shown  in  Fig.  16 
(15  UTC  in  Fig.  16a,  16  UTC  in  Fig.  16b,  and  17  UTC  in  Fig.  16c).  The  modeled 
convection  in  southeastern  Kansas  evolves  closer  to  the  observed  linear  structure 
(Fig.  16d)  at  later  forecast  start  times,  while  the  modeled  convection  just  south  in 
northeastern  Oklahoma  strengthens  contrary  to  the  observed  gap  in  the  line  here. 
Encouragingly,  the  convection  just  north  of  Moore  strengthens  considerably  with 
shorter  lead  times,  so  that  in  OBSGRID3O20,  greater  than  70  dBZ  is  simulated 
north  of  Moore  (Fig.  16c).  Although  the  strongest  modeled  convection  is 
somewhat  north  of  the  observed  storm,  this  solution  is  much  closer  to  the 
observed  radar  reflectivity  than  earlier  cycles.  By  the  OBSGRID3O20  cycle 
(Fig.  16c),  a  broken  line  of  precipitation  in  central  Kansas  forms  near  the  axis  of 
the  observed  line  in  this  region  (Fig.  16d),  but  the  modeled  line  is  not  as 
continuous  and  not  as  wide  as  the  observations.  In  addition,  behind  this  line 
OBSGRID3O20  erroneously  shows  a  strong  line  of  convection  stretching  into  far 
northern  Texas.  0BSGRID30i9  and  OBSGRID3O20  also  show  stronger 
convection  in  the  southern  quarter  of  the  domain  than  observed.  Also,  the 
erroneous  area  of  convection  in  western  Oklahoma  simulated  in  many  of  the 
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experiments  becomes  much  more  prominent  with  shorter  lead  time.  In  general, 
0BSGRID30  appears  to  simulate  overly  strong  convection  at  cycles  with  shorter 
lead  times  (e.g.,  OBSGRID3O20)  but  better  simulates  the  convection  in 
southeastern  Kansas  and  in  the  Moore,  Oklahoma,  region. 


Fig.  16  WRF-ARW  lowest  model  level  radar  reflectivity  at  2015  UTC  on  20  May  2013  for 
a)  0BSGRID30i8,  b)  0BSGRID30i9,  c)  OBSGRID3Cho,  and  d)  observed  base  reflectivity 
composite.  The  location  of  the  Newcastle-Moore  tornado  at  2015  UTC  is  indicated  by  the 
black  “X”  within  the  white  square. 

The  full  hybrid  experiment  (VLAPS3AO)  includes  both  the  vLAPS  analysis  in 
initial  conditions  as  used  by  VLAPS3,  the  observation  nudging  used  in 
0BSGRID30  and  VLAPS30,  and  the  analysis  nudging  toward  vLAPS  analyses 
during  the  first  3  h  used  in  VLAPS3A.  VLAPS3AO  simulates  the  northeastern 
part  of  the  main  line  of  convection  (in  southeastern  Kansas)  for  VLAPS3AOi8 
(Fig.  17a),  VLAPS3AOi9  (Fig.  17b),  and  VLAPS3A02o  (Fig.  17c),  albeit  with 
more  scattered  convection  along  the  axis  of  the  line  than  observed.  VLAPS3AOi9 
and  VLAPS3AO20  forecast  less  convection  in  the  portion  of  the  main  convective 
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line  in  far  northern  Oklahoma  than  VLAPS3AOis  and  more  consistent  with  the 
gap  seen  in  the  line  in  the  2015  UTC  radar  reflectivity.  Overprediction  of 
convection  in  the  southern  quarter  of  the  domain  in  VLAPS3AOi8  and 
VLAPS3AOi9  is  improved  in  VLAPS3AO20.  Overall,  the  VLAPS3AOi9  cycle  in 
particular  overpredicts  the  coverage  of  convection.  While  VLAPS3AOis  predicts 
west-to-east  convection  just  north  of  the  Newcastle-Moore  tornado, 
VLAPS3AOi9  and  VLAPS3AO20  both  predict  strong  convection  (~70  dBZ)  just 
north  of  the  tornado  as  part  of  a  southwest-northeast  oriented  element.  These 
cycles  in  the  experiment  using  vLAPS  analysis  for  initial  conditions  without  any 
nudging  (VLAPS3i9,  Fig.  14b;  VLAPS320,  Fig.  14c)  did  not  produce  convection 
close  to  the  Newcastle-Moore  tornado  as  strong  as  that  in  YLAPS3 AO. 
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Fig.  17  WRF-ARW  lowest  model  level  radar  reflectivity  at  2015  UTC  on  20  May  2013  for 
a)  VLAPS3AOis,  b)  VLAPS3AOi9,  c)  VLAPS3AO20,  and  d)  observed  base  reflectivity 
composite.  The  location  of  the  Newcastle-Moore  tornado  at  2015  UTC  is  indicated  by  the 
black  “X”  within  the  white  square. 
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6.2  Ability  of  Hybrid  to  Simulate  Strong  Moist  Convection 


While  the  previous  section  focused  on  comparing  different  experiments  and 
cycles,  here  we  demonstrate  the  capability  of  the  full  hybrid  experiment 
(VLAPS3AO)  to  simulate  strong  convection.  The  goal  here  is  not  to  demonstrate 
that  VLAPS3AO  simulates  the  supercell  that  spawned  the  Newcastle-Moore 
tornado  consistently  with  observations,  but  rather  to  show  that  while  using  the  full 
hybrid  data  assimilation  technique,  the  model  can  simulate  a  strong  storm  with 
features  consistent  with  those  expected  in  strong  moist  convection. 

The  model-simulated  radar  reflectivity  at  the  lowest  model  level  is  shown  for 
VLAPS3AOi8  at  2015  UTC  in  Fig.  18  along  with  the  model-simulated  wind  at  the 
lowest  prognostic  level  (~13  m  AGL).  Ahead  of  the  main  line  of  convection, 
winds  are  from  the  south  while  behind  the  line,  there  is  a  lot  of  variation  in  the 
wind  directions.  In  order  for  the  wind  barbs  to  be  legible,  only  winds  from  a  small 
portion  of  the  model  grid  points  are  plotted  (the  model  has  800  grid  points  in  each 
direction).  We  focus  on  a  storm  southeast  of  the  actual  location  of  the  Newcastle- 
Moore  tornado  (small  white  square  with  black  X) — namely,  looking  at  the  area 
covered  by  the  white  box  with  horizontal  plots  and  at  the  west-east  line  (or 
portions  thereof)  with  cross  sections.  The  zoomed-in  view  of  the  lowest  model 
level  reflectivity  and  winds  (Fig.  19)  shows  a  core  of  approximately  70  dBZ 
reflectivity  within  a  much  larger  area  of  40+  dBZ  reflectivity.  A  much  more 
detailed  view  of  the  surface  winds  is  available  in  this  zoomed-in  image  (Fig.  19) 
compared  to  the  whole  domain  view  (Fig.  18).  Here,  south  winds  are  seen  at  the 
surface  outside  of  the  storm  (southeast  comer),  with  winds  flowing  strongly  away 
from  the  area  with  highest  reflectivity  except  on  its  western  side  (the  wind  barbs 
indicate  winds  of  25+  ms'1  on  the  east  side  of  this  area).  An  area  with  slightly 
lower  model-simulated  reflectivities  (peaking  in  the  mid-60s  dBZ)  is  predicted 
just  west  of  the  area  with  maximum  reflectivities. 


Approved  for  public  release;  distribution  is  unlimited. 


49 


102°W  100°W  98°W  96°W 


Fig.  18  WRF-ARW  lowest  model  level  radar  reflectivity  at  2015  UTC  on  20  May  2013  for 
VLAPS3AOi8.  The  location  of  the  Newcastle-Moore  tornado  at  2015  UTC  is  indicated  by  the 
black  “X”  within  the  small  white  square.  The  large  white  square  indicates  the  area  shown  in 
Fig.  19,  and  the  white  horizontal  line  indicates  the  location  of  the  cross  section  shown  in 
Fig.  20.  Each  half  barb  is  2.5  ms'1  and  each  full  barb  is  5.0  ms-1. 
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Fig.  19  WRF-ARW  lowest  model  level  radar  reflectivity  at  2015  UTC  on  20  May  2013  for 
VLAPS3AOi8  zoomed  into  the  area  denoted  by  the  white  box  in  Fig.  18.  Each  half  barb  is 
2.5  ms-1  and  each  full  barb  is  5.0  ms-1. 

A  vertical  cross  section  of  the  storm  (Fig.  20)  shows  2  narrow  cores  of  high 
reflectivity  (~70  dBZ)  and  some  graupel.  The  stronger  and  farther  eastern  core  has 
strong  vertical  motion  just  east  of  it;  the  area  of  graupel  includes  the  area  above 
this  updraft.  An  anvil  is  also  predicted  stretching  well  east  of  the  horizontal  extent 
of  the  storm  at  ground  level.  It  is  difficult  to  see  the  full  model  structure  of  the 
main  storm  in  this  cross  section  since  for  clarity  wind  is  plotted  at  many  fewer 
grid  points  than  are  present  in  the  model  (horizontally  wind  vectors  are  plotted 
every  ~5  grid  points). 
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Fig.  20  Cross  section  of  VLAPS3AOis  along  the  white  horizontal  line  in  Fig.  18  showing 
winds  in  the  plane  (i.e.,  the  model  u-  and  w-component  winds  since  this  cross  section  is  along 
the  model  x-dimension),  model  simulated  radar  reflectivity  (shaded  with  a  thick  white  line 
indicating  the  0  dBZ  isoline  and  a  thin  white  line  indicating  the  -25  dBZ  isoline),  and 
graupel  mixing  ratio  (black  contours  every  3  g  kg-1).  The  thick  black  line  at  the  bottom  of 
the  plot  indicates  the  portion  of  the  cross  section  included  in  Figs.  21  and  22.  Each  half  barb 
is  2.5  ms-1,  each  full  barb  is  5.0  ms1,  and  each  flag  is  25.0  ms1. 

The  model-predicted  structure  of  this  storm  is  seen  more  clearly  in  the  zoomed-in 
cross  section  in  Fig.  21.  Near  the  surface,  just  east  of  the  highest  reflectivity 
values  (i.e.,  approximately  model  grid  cell  32  in  terms  of  the  x-axis  label),  the 
model  predicts  strong  west  winds  (=28  ms-1),  consistent  with  an  outflow.  About 
10  km  farther  east  (e.g.,  model  grid  cell  42),  the  flow  in  the  transect  just  above  the 
surface  reverses  from  west  to  east,  with  the  east  wind  serving  as  the  inflow  to  a 
strong  updraft  between  model  grid  cells  38 — 42.  The  updraft  includes  a  region  with 
weaker  echoes  than  found  to  the  west  and  to  the  east,  but  the  magnitude  of  the 
echoes  is  still  notable.  By  7,000  m  AGL  in  the  updraft,  the  reflectivity  increases 
to  approximately  65  dBZ  with  an  area  of  graupel  and  graupel  mixing  ratios 
exceeding  6  g  kg  1  around  8,000  m  AGL.  Just  to  the  west  of  the  column 
containing  the  strong  updraft  with  graupel  in  its  upper  portion  is  a  column  with 
strong  reflectivity  reaching  down  to  the  surface.  For  grid  cells  28-32, 
approximately  70  dBZ  reflectivities  reach  the  surface  and  3g  kg  1  graupel  mixing 
ratios  reach  as  low  as  2,250  m  AGL  (graupel  mixing  ratios  of  >9  g  kg"1  are 
forecasted  around  7,000  m  AGL  near  grid  cell  32)  .  Reflectivities  greater  than  or 
equivalent  to  60  dBZ  cover  most  of  the  column  between  0  and  7,000  m  AGL  for 
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grid  cells  26-32.  The  main  updraft  in  the  storm  is  strong  (Fig.  22),  peaking  at 
greater  than  35  ms-1  around  6,700  m  AGL  around  grid  cell  41.  There  are  also 
strong  updrafts  predicted  at  the  top  of  the  main  column  with  high  reflectivities; 
updrafts  peak  at  greater  than  35  ms  1  around  9,000  m  AGL  around  grid  cell  34. 
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Fig.  21  Cross  section  of  VLAPS3AOis  along  the  portion  of  the  cross  section  in  Fig.  20 
indicated  by  the  thick  black  line  along  the  bottom  of  that  figure.  This  cross  section  shows 
winds  in  the  plane  (i.e.,  the  model  u-  and  w-component  winds  since  this  cross  section  is  along 
the  model  x-dimension),  model-simulated  radar  reflectivity  (shaded  with  a  thick  white  line 
indicating  the  0  dBZ  isoline  and  a  thin  white  line  indicating  the  -25  dBZ  isoline),  and 
graupel  mixing  ratio  (black  contours  every  3  g  kg-1).  Each  half  barb  is  2.5  ms-1,  each  full 
barb  is  5.0  ms-1,  and  each  flag  is  25.0  ms-1. 
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Fig.  22  Cross  section  of  VLAPS3AOis  along  the  portion  of  the  cross  section  in  Fig.  20 
indicated  by  the  thick  black  line  along  the  bottom  of  that  figure.  This  cross  section  shows 
model-simulated  radar  reflectivity  (shaded  with  a  thick  white  line  indicating  the  0  dBZ 
isoline  and  a  thin  white  line  indicating  the  -25  dBZ  isoline),  and  vertical  motion  (black 
contours  every  5  ms-1  with  contours  representing  downward  motion  dashed).  Vertical 
motion  ranges  from  -16  to  +38  ms-1  in  this  cross  section. 

VLAPS3AOi8  is  able  to  forecast  strong  moist  convection,  with  elements  consistent 
with  those  expected  to  be  present  in  strong  convection.  The  reflectivities,  horizontal 
wind  speed,  updrafts,  and  graupel  predicted  by  the  model  are  consistent  elements 
expected  in  a  strong  thunderstorm.  This  example  demonstrates  the  capability  of  this 
WRF  configuration  to  simulate  strong  moist  convection. 

7.  Summary,  Discussion,  and  Conclusions 


A  preliminary  analysis  of  WRF  simulations  using  vLAPS  analysis  in  combination 
with  observation  and  analysis  nudging  data  assimilation  suggests  potential  value 
to  this  new  combination  data  assimilation  method  when  evaluated  for  a  single 
case.  The  vLAPS  nudging  methodology  combines  2  data  assimilation  techniques 
both  with  limited  computational  requirements.  The  methodology  leverages  the 
ability  of  vLAPS  to  create  analyses  containing  in  situ  and  remotely  sensed  data 
with  the  ability  of  nudging  to  assimilate  continuously  and  gradually.  While 
observation  nudging  can  assimilate  observations  of  fields  that  are  prognostic 
fields  within  the  model,  it  cannot  directly  assimilate  observations  of  fields  such  as 
radar  reflectivity  that  are  not  model  prognostic  variables;  vLAPS  can  assimilate 
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fields  whether  or  not  they  are  model  prognostic  variables.  Although  vLAPS 
accounts  for  temporal  variability  among  the  observations,  the  resultant  analyses 
are  normally  directly  inserted  as  the  model  initial  conditions;  analysis  nudging 
allows  these  analyses  to  gradually  affect  the  model  solution  as  well  as  allowing 
the  assimilation  of  a  series  of  vLAPS  analyses  in  a  single  model  integration.  The 
gradual  assimilation  of  a  series  of  vLAPS  analyses  should  allow  the  model  initial 
conditions  at  the  end  of  the  assimilation  to  be  more  consistent  with  the  model 
dynamics  and  physics  and  thus  more  likely  to  be  retained  by  the  model. 

In  this  study,  the  model  experiments  differed  in  the  source  of  the  initial 
conditions,  whether  analysis  nudging,  observation  nudging,  or  neither  were 
applied,  and  the  amount  of  model  integration  time  before  the  time  designated  as 
the  0-h  forecast  (i.e.,  the  length  of  the  model  spin-up).  For  the  initial  conditions,  a 
single  cycle  of  the  3-km  HRRR  was  either  used  directly,  or  used  as  the  first-guess 
field  for  either  a  vLAPS  analysis  (variational  for  some  fields)  or  an  analysis 
performed  by  Obsgrid  (modified  multiscan  Cressman  analysis).  Experiments 
using  analysis  nudging  nudged  toward  vLAPS  analyses,  where  experiments  using 
observation  nudging  nudged  toward  various  MADIS  observation  types. 
Additionally,  multiple  cycles  of  each  experiment  were  performed,  with  a  1-h 
cycle  interval  and  each  cycle  having  no  influence  on  the  subsequent  cycle.  Given 
the  9  experiments  and  the  6  cycles  per  experiment,  there  are  a  large  number  of 
model  integrations  to  evaluate,  even  accounting  for  the  model  simulations  that 
served  as  one  cycle  in  one  experiment  and  as  another  cycle  in  another  experiment. 
The  large  amount  of  model  output  to  evaluate  requires  an  objective  methodology, 
but  here  we  performed  a  preliminary  subjective  analysis  by  focusing  on 
comparing  model-predicted  reflectivity  to  observed  reflectivity  at  2015  UTC.  The 
primary  goal  is  to  evaluate  whether  a  hybrid  data  assimilation  approach 
combining  analysis  nudging  toward  vLAPS  analysis  and  observation  nudging 
performs  better  than  either  vLAPS  or  nudging  alone. 

The  full  hybrid  experiment  (VLAPS3AO)  showed  promise  in  its  radar  reflectivity 
predictions  as  compared  to  other  experiments.  There  were  indications  that 
providing  a  longer  time  period  for  the  model  to  spin  up  can  be  beneficial.  Overall, 
it  was  difficult  to  determine  from  this  limited  subjective  comparison  which 
experiments  performed  best.  It  should  also  be  noted  that  by  looking  at  a  single 
output  time,  we  did  not  consider  the  potential  that  experiments  might  be 
producing  convection  that  appears  to  be  a  poor  forecast  at  the  time  evaluated 
(here  2015  UTC)  due  to  the  model  forecasting  a  convective  element  too  early  or 
too  late.  The  general  consistency  among  the  experiments  in  the  overall  forecast 
for  moist  convection  (e.g.,  the  strong  south-west  to  north-east  line  of 
convection),  combined  with  the  significant  variations  among  the  forecasts  on  the 
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details  of  the  convection,  appears  consistent  with  the  work  of  Zhang  et  al.  (2015, 
2016)  examining  predictability  for  this  case.  A  more  quantitative  and  thorough 
analysis  of  our  WRF  simulations  is  needed  to  determine  the  potential  value  of  the 
different  assimilation  techniques,  in  light  of  the  limits  of  predictability  suggested 
by  Zhang  et  al.  (2015,  2016). 

We  also  demonstrated  the  capability  of  VLAPS3AO  to  simulate  strong  moist 
convection.  A  storm  predicted  by  VLAPS3AO  was  discussed  that  included 
various  elements  associated  with  strong  moist  convection.  WRF  predicted 
70+  dBZ  reflectivities,  35+  ms'1  updrafts,  and  25+  ms"1  outflow  winds  at  the 
surface,  along  with  predicting  graupel  in  the  storm. 

A  more  objective  and  complete  analysis  is  required  to  better  characterize  the 
relative  value  of  the  different  assimilation  techniques  being  evaluated  in  this 
study.  Evaluating  all  cycles  of  all  experiments  at  all  forecast  times  objectively 
will  allow  for  more  definite  conclusions  to  be  drawn.  Even  then,  it  must  be 
remembered  that  this  is  a  single  case,  and  future  work  should  choose  those 
experiments  showing  the  best  potential  and  evaluate  them  for  further  cases, 
including  cases  lacking  the  strong  convection  that  characterized  this  case.  Also, 
the  evaluation  of  reflectivity  fields  is  particularly  difficult,  given  the  strong 
gradients  in  the  field,  especially  in  convective  precipitation.  Object-based 
methods  can  provide  improved  insight  during  the  evaluation  of  radar-reflectivity 
forecasts  but  can  be  challenging  and  time  consuming  to  apply  well. 

In  addition  to  evaluating  additional  case  days,  other  areas  of  future  research 
would  also  be  beneficial.  The  analysis  nudging  toward  vLAPS  analyses  is 
currently  limited  to  the  fields  WRF  has  the  capability  to  nudge  toward  (potential 
temperature,  water  vapor  mixing  ratio,  and  winds).  Future  work  could  explore 
analysis  nudging  toward  other  fields  in  the  vLAPS  analyses  (e.g.,  cloud  water 
mixing  ratio)  to  more  fully  utilize  the  vLAPS  analyses.  Other  sources  of  analyses 
to  nudge  toward  could  also  be  explored. 

This  report  describes  the  combination  of  nudging  and  vLAPS  data  assimilation 
techniques  and  provides  preliminary  results  from  the  application  of  this  technique 
to  a  single  case  day.  The  combination  technique  shows  promise  in  that  it  can 
assimilate  observations  of  both  prognostic  WRF  fields  (e.g.,  water  vapor, 
potential  temperature,  and  winds)  and  other  fields  (e.g.,  radar  reflectivity)  and  can 
do  so  in  a  framework  that  gradually  applies  the  effects  of  these  observations  to  the 
model  solution  which,  it  is  hoped,  will  allow  the  model  to  better  retain  the 
analyses  than  if  they  were  inserted  into  the  model  at  a  single  time  step.  Further 
analysis  of  this  case  is  needed  to  better  understand  the  value  of  the  combination, 
and  further  case  studies  are  needed  to  confirm  results  from  the  current  study. 
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List  of  Symbols,  Abbreviations,  and  Acronyms 


1-D 

1 -dimensional 

3-D 

3 -dimensional 

4-D 

4-dimensional 

3DVAR 

3 -dimensional  variational 

4DVAR 

4-dimensional  variational 

ABL 

atmospheric  boundary  layer 

ACARS 

Aircraft  Communications  Addressing  and  Reporting  System 

AGL 

above  ground  level 

ARL 

US  Army  Research  Laboratory 

ARPS 

Advanced  Regional  Prediction  System 

EnKF 

ensemble  Kalman  filter 

FGAT 

first  guess  at  appropriate  time 

GRIB 

Gridded  Binary 

GSI 

Gridpoint  Statistical  Interpolation 

HNEnKF 

hybrid  nudging-ensemble  Kalman  filter 

HRRR 

High  Resolution  Rapid  Refresh 

LAPS 

Local  Analysis  and  Prediction  System 

MADIS 

Meteorological  Assimilation  Data  Ingest  System 

MYJ 

Mellor- Y  amada- Janj  ic 

NCAR 

National  Center  for  Atmospheric  Research 

NetCDF 

network  common  data  format 

NWP 

numerical  weather  prediction 

TKE 

turbulent  kinetic  energy 

UTC 

Coordinated  Universal  Time 

vLAPS 

variational  Local  Analysis  and  Prediction  System 

WPS 

WRF  Preprocessing  System 
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WRF 

WRF-ARW 

WSR-88D 


Weather  Research  and  Forecasting  model 

Advanced  Research  version  of  the  Weather  Research  and 
Forecasting  model 

Weather  Surveillance  Radar,  1988,  Doppler 
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